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ABSTRACT 


This thesis consists of an interactive program that enables the student to study the 
orbital motion of satellites around the earth. The student can investigate the shape of 
a Variety of orbits by varving the initial position and velocity of the satellite. or by sup- 
plving select orbital parameters 1.e. initial orbital radius, eccentricitv, and inclination. 
Satellite maneuvers can also be studied, like transfer orbits and inclination changes, by 
command velocity changes at any location in the orbit. Also the effects of the perturb- 
ing forces due to the oblateness of the earth. drag for low earth orbits, and gravitational 
attraction from the sun and moon can be investigated. The orbits are displaved in either 
the perifocal coordinate system around a model of the earth. or the ground track can 
be displaved on a map of the world. Orbital data is displaved below the orbital plot. 
The display is enabled by the use of display integrated software svstem and plotting 


КЫ злее (DISSPLA) subroutines. 
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1. INTRODUCTION 


A visual aid for students new to orbital mechanics is required to comprehend fully 
the dynamics of orbital motion. This program is an interactive time step simulation 
program that calculates and plots either unperturbed or perturbed elliptical orbits. The 
program interacts with the student in developing the iniual orbit. Also the program 
emeoles the student with the ability to change the velocity of the satellite at a specific 
location in the orbit. This feature will permit the student to investigate the effects of 
commanded velocity changes as in perigee kicks, apogee kicks and inclination changes. 
The user can also modifv the iniual position and velocity of the satellite at the com- 


pletion of anv orbit. 


The student is given an opportunity to investigate the effects of perturbing forces 
on the satellites orbit by choosing to have the program calculate the orbit with or with- 
ШО Огоо forces. The variation of parameters method, as seen in [Ref. 1: рр. 
396-307]. 1s used in calculating the perturbing orbit. The perturbing forces taken into 
consideration are the following: 


1. the oblateness of the earth 


Коле [ог [ом carin orbits 
anatona] force of the moon 
4. сгаупапопа force of the sun 


In order to review fullv the operation of the program (included in appendix A) and 
to uncover anv problems or limitations that plagued the programming. the program has 


been divided up as follows: 


l. program design 


tly 


unperturbed orbit 


(5,2 


perturbed orbit 


4. velocity changes 


сл 


graphicai plots 


Ihe programming approach and equations used in each of the above sections wil! be 


оси а there respective ciapters. A review of the coordinate systems used and their 


transiormauons between them are meluded in appendix B. Since all the cquatons teem 
in the calculation of the orbital elements are from reference |, they will not be reviewes 
in each chapter but will be included in appendix C fora quick zeference. Laquationsmeen 
other sources. will be referenced in Оет осе e 

Examples of perturbed and unperturbed orbital plots for a variety of initial orbital 
paramicters are included in appendix D. Included are plots of low earth orbits, transfer 


orbits and geosvnchronous orbits. 


Ill. PROGRAM DESIGN 


In designing this program an attempt was made to make it not onlv as user friendly 
as possible, but also to make the program as simple as possible to understand. To 
achieve these goals, the program would have to be written in a logical manner, in a 
computer language that 1s easy to follow, the program would have to run on terminals 
readilv available to students (at the Naval Postgraduate School (NPS)), and the program 
would have to be easily used by students with a minimum amount of computer or orbital 


mechanics knowledge. 


FORTRAN was chosen as the programming language since it is a. wildlv used sci- 
entific language and it allows for very structured programming. By programming in a 
structured format. the program can be expanded in the future with a minimum amount 
of time required to understand the programming code. FORTRAN also allows for 
double precision numbers to be used in the calculation of the orbit. This ts critical when 
round off error in single precision could be greater then the actual change that one 1s 
trving to model. The equations in the descriptions of the program might not exactly 
match the equations in the listings because of special programmung techniques which 
must be included in most computer programs to handle such problems as “division by 


zero . 


The displav integrated software svstem and plotting language (DISSPLA) package 
available on the mainframe computer at NPS was used to enable a variety of graphical 
displays with a minimum amount of programming. DISSPLA has a set of subroutines 
that the programmer calls to display data contained in arrays. This requirement forces 
the program to load arravs with the satellites position in order for it to be plotted. The 
TEC618 computer terminal and associative plotter was used for ease of gaining hard 
gop plots of the orbits and the diversity of locations that are available here at NPS. 
In order to run a program in DISSPLA the user must first define storage space of 1300k 
and designate temporary disk space, and then call DISSPLA with the program name. 


This 1s accomplished with the following commands: 


ENDCLINESIORAGE 1500K 


UNE 
>. ТРЕЕ Dis 
вЫ 


To make the program user friendly, the user is prompted for inputs via the keyboard. 
The entry is usually a number. X ves or no response can be entered bv tvping НЕ 
“№”. Та most cases the program does a check to see if the input 1s appropriate. In order 
to make it as easv as possible for the student to get the desired orbit displaved, the 
program requires onlv the initial position and velocity of the satellite. The initial posi- 
tion and velocity of the satellite 1s supplied by the user in one of two ways. The user can 
input the position and velocity of the satellite. using the perifocal coordinate system 
(ОК). ог the user can let the program place the satellite on the 1° axis of the IJK system 
at the radius of perigee (RP) distance supplied by the user. This latter choice gives the 
initial location of the satellite, but to get the velocitv the program will prompt the user 
for one of the following: 


l. -the actual veloci асе ES STEM 


2. the eccentricitv (e) of the orbit. In which case the velocity 15 calculated поши 
following equations: 


[Р 


1-- е 





C= = semi-major aX1s 
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ENR=- = energy mass 
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о, 
NI mass of earth 
G = Lniversal gravitational constant 
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; ц 
. = Й У - и -- pep 


the radius of apogee (RA). The velocity is calculated by first calculaung them 
centricitv (e) from the following: 


_ ВА- ВР 
Е 


With the eccentricity the same equations used above are used to calculate the ve- 
loe 


(22 


е 


In order to give the velocity a direction the inclination (1) of the orbit 1s required from 


the user. The following equations are used to calculate the velocity vector: 


уу = 0.0 


I 


v = Y COS(I) 


v, 25 v Sint) 


The program will check to ensure that the orbital eccentricity is less than 1.0, if it is not 
then the program will reject the inputs. After the initial input are accepted, the program 
will do calculations for the six orbital elements required to describe the size. shape and 
Orientation of the orbit, and to pinpoint the position of the satellite along the orbit at a 
particular ume. This classical set of six orbital elements are as follows: 

|. a. semi-major axis. 

с eccentricity. 


1, inclination. 


„ә 


62, longitude of the ascending node. 


Е 


«), argument of perigee passage. 


m 


емле о? регсее passage. 


The program actually calculates more orbital elements than the six classical elements 
required to plot the orbit, this 1s done in an effort to make the program as robust as 


possible. This will add in the ability to expand the program in the future. 


ete Satellite 1s not tninally at the perigee point then the satellite must first be 
stepped around to the perigee point. The program then enters a loop that calculates the 
orbit from the perigee point through one complete orbit around the earth and back to 
the perigee point. The orbit ts calculated in steps of 2 times pi divided bv an integer, Le., 
2 umes pi divided bv 50. This step size was used to ensure a smooth orbit for display 
purposes and also to get within adequate distance to the perigee point or other location 
for a velocity change. After the loop 1s completed, the program will offer the user a 
choice of the following plots to check the orbit: 


1. perifocal 


2. groundtrack 


The program then goes into a loop offering the user the following choices until the user 


Сес то end the prosran 


l. plot another view of the same orbit. 
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IHI. UNPERTURBED ORBIT 


The subroutines that calculate the unperturbed orbit are the most widely used sub- 


routines in the entire program. These subroutines are called to step the satellite around 


to the perigee point from the user supplied initial position and velocity, to calculate the 


next unperturbed orbit, and for anv velocity change. No matter wnich of these sources 


Supply the initial position and velocity the program calculates the unperturbed orbit in 


the same manner. The only difference 1s where in the orbit the satellite is initially when 


these subroutines are called. Before the unperturbed subroutines are called, the orbital 


elements are calculated. 


The unperturbed subroutines are called bv a single subroutine L NPRET’ which has 


the following basic algorithm: 


Г. 


Increment time by the time step size (DT). The time step was chosen as the period 
mica by hfiv to give a smooth plot. but more importantly to ensure that the 
satellite 1s within an acceptable distance from a specific location for a velocity 
enamce. Nhe angular error caused by the step size can be as much as PI 50 from 
the desired point for a circular orbit and will increase for more eccentric orbits. 
This error becomes a factor when the user 1s making velocity changes. and therefore 
it Will be covered in that chapter in further detail. 


аге the new elements. The calculation of the new elements ts the heart of this 
algorithm. The size, shape and orientation of the orbit remains unchanged. What 
is required is the position of the satellite along the orbit as a function of time. The 
problem becomes a matter to solve ^ = Kepler problem -predicung the future po- 
sition and velocity of an orbiting object as a function of some known iniual posi- 
tion and velocitv and the time of fii aif Дер ои using these 
principles will follow: 


a. A ume step (DT) is added to the time of flight( TF). time of flight is the elapsed 
time since the satellite passed the perigee point. 


ша iP ОГ 


b. The new mean anomaly (MA) 1s calculated from the new time of flight, and the 
mean moton (MM). 


МА = ММ х ТЕ 


c. With the new mean anomaly the new eccentric anomaly (EA) is calculated. 
селе solution (0 the Kepler problem (MA - E4 —exsm(Ed)) ts 
uanscendental, an iterative solution based on the Newton method of root find- 
Inga is. (сес. The root in question is a solution to the equation 
(Ma — E.f +e x sin( EA) =0). This algorithm takes the form of [Ref. 1: p. 222]: 


ПИ лм гк Е) 


{ 2 


E к i i-r 
инь: = 21 (l—e x cost£4,)) 


Where this equation 1s applied initially to EA, = МАЯ апа then reapplied 
until the difference between MA and M4, becomes small enough to be ig- 
nored. 


d. The new true anomali и је о ва Јаке on 


cos (e — cos( EA )) 
ше е cos(EA) — 1 


(+2 


Calculate the new position and velocity. The position and velocity are calculated 
in the perifocal coordinate svstem (PQW). The PQW svstem uses the orbit as its 
fundamental plane and therefore requires only two coordinate to specify thet amii 
ше 5 position and velocity. The z, coordinate is bv definition always equal to zero. 
The position of the satellite is calculated as: 


xX, =rcos vy 
== р] ү 


= () 


“= 


The velocity of tie аге ВО as: 


l ! 
b Se E ( =з ге) 
Шан 
v = (ей со ту) 


M 


y=) 


Ф 


4. Store position and elements in arrays for plotting. In order for the program to plot 
the orbit the radius, true anomaly. inclination, and argument of perigee must be 
stored in arrays. The use of these arrays to plot the orbit will be explainedaim 
chapter 6. 


5. The process 1s repeated until the satellite 1s at the perigee point and themes 
anomaly 1$ two pi. 

The procedure used to calculate the unperturbed orbit leave verv little to be modified 
by a programmer. The only choices that had to be made concerned step size, how to tell 
the UNPRET subroutine that tlie pengee point пай been reached. ane@™a заме она 
ceptable error for newtons method. For the unperturbed orbit, the step size just had to 
be small enough to produce a smooth plot of the orbit. Two indicators for perigee were 
used, one was that the true anomaly was ereater than 6.21 radians (two pi equals 6.25 
radians) and the time from the previous perrgee point will be greater then the period. 


The two indicators were logically ‘and’ together to ensure the perigee point was reached: 


The disparity between two pr and 6.21 radians 1s due to the error produced bx the sat- 
mere wot becinning the orbit at exactly the perigee point and the step size used go 
пошта the orbit. The acceptabie size of error for ne«tons method was set at. 1.01I:-10, 
because for an unperturbed orbit this would be the major contributor to anv error in the 
ЕШ апа the magnitude of this error would be acceptable. However: in a perturbed 
orbit there are other factors contributing to determining the acceptable error, and these 


will be discussed in the next chapter. 


IV. PERTURBED ORBIT 


The perturbed orbit uses the same basic routines as the unperturbed orbit 1п згер- 


ping 


the satellite around the earth with one major difference, the perturbing forces 


produce a time rate of change of the orbital elements that must be applied at each time 


step. 


The variation of parameters method is used to determine this influence C Emig 


perturbing forces on the orbital elements. The analvsis is simplified bv using the orbital 


Coor 


dinate svstem ‘RSW’, as explained in appendix B. The basic algorithm 1s as follows 


[Ref. 1: p. 407]: 


I Atr =n calculate six OrDita CEDE 
2. Compute the perturbing forces and transiorm 1с асг = г to the RSW' SYSTEMB 
5. Compute the time rate-of-change of the elements. 
4. Calculate the change of elements for one time step, and add the changes to the old 
хайшез ат еасй ер то сег тие пе Cue ES. 
5. From the new values of the orbital elements. calculate a position and velocity. 
6. Goto the step and repeat until the inal mme тел По 
The steps in the algorithm will be explained in the following sections: 
A. ORBITAL ELEMENTS 
The standard orbital elements ences) апл ое 
a = semi-major axis 
e = eccentrici 
1 = inclination 
Q = longitude of ascending node 
© = argument of perigee 
T = time of perigee passage 
(M, = mean anomalv at epoch » M — n(r— t))). The elements are calculated only 


at the beginning of the orbit from the initial position and velocity vectors. The elements 


are 


then changed continuouslv throughout the orbit bv adding the changes due to the 


perturbing forces. For the perturbed orbit, the satellite will always begin at the perigee 


point. This is done so one complete orbit 15 from perigee point to perigee point. 
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ШИСОМРЕТЕ РЕКТСВКВГУС FORCES 

The variation of parameters method requires that the perturbing forces be calculated at 
each step in the orbit. In order to do this a model of each perturbing force must Бе 
developed. The following perturbing forces where used in calculating the total perturb- 
ing force effecting the satellite: 


1. oblateness of the earth 


ty 


atmospheric drag 


G- 


gravitational attraction of the sun 


4. gravitational attraction of the moon 


The magnitudes of these forces have an enormous range of values and are dependent 
on the distance the satellite is from the perturbing body. Figure 1 on page 12 shows a 
graphical representation of the magnitude of the perturbing forces in a log-log plot of 
Pemunbinegeiorces per unit mass [Ref. 2: p. [V-61}. The model of each of these forces 
follows: 

1. NON-SPHERICAL EARTH 

The earth is not perfectly spherical, but bulges around the equator. The polar 

amsesuatoral diameters are 12713.0 km and 12/56.5 Km. respectively. The oblateness 
results in a perturbing force per unit mass with these components in the 'RSW' coordi- 
nate svstem [Ref. 3: p. 81]: 


F, 2 ——Fc— (1-3 sin'(i) sin (9) 


қ 
С 3 


ШЕ (sin (i) sinko) cosig)) 
р 
( Е ЦИ) rr Ар. 
p — — (sin(i) cos(?) sin(ug)) 


The variable and constants of these equations are defined below: 
Ши маша! ез: 


а. и, = the argument of latitude and is equal to the true anomaly v, plus the ar- 
cument of perigee c. 


ио = 7 = 0 


b. r = the radius from the center of the earth to the satellite. 


1] 
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Figure 1. Comparison of perturbation magnitudes. 
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EN onstarnts: 


а. u = the gravitational parameter of the earth, 
3 
В (kin) 
i SSO CE 
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b. J; = the second harmonic of oblateness coefficient, determined by experimental 
observations, 


J, = 1.0823E — 3 


~ 


с. г = Ще mean radius of the earth, 


e 


r, = 6.3782 E3Km 


2. ATMOSPHERIC DRAG 
The formulation of atmospheric drag equations are plagued with uncertainties 
of atmospheric fluctuations, frontal areas of orbiting object (if not constant), the drag 
Seer ieient, ald other parameters. A fairl simple formulation will be given here. Drag, 
by definition, will be opposite to the velocity of the vehicle relative to the atmosphere. 


Thus. the perturbing force 1s 


Е = (5). Ср. лк. ЕХ. х. 


~ 


Ши оси хестог 15 1 е ШК зузсет 50 ше гезшипо force is also in the IJK’ sys- 
tem. Therefore a transformation to the 'RSW' svstem is required. 
The variables and constants of this equation are defined below: 
]. Variables: 
me = speed of vehicle. 


b. CD = the dimensionless drag coefficient. The drag coefficient CD has a value 
between | and 2. It takes a value near 1 when the mean free path of the at- 
mospheric molecules 1s small compared with the satellite size, and takes a value 
close to 2 when the mean free path 1s large compared with the size of the satel- 
hte. The drag coefficient will be modeled with CD = 2 when the satellites alti- 
tude is greater than 550km and equal to | otherwise. [Ref. 4: p. 295] 


c. DEN = atmospheric density at the vehicle’s altitude. The density 1s spherically 
symmetric, and will be modeled using exponential steps using the parameters in 
Table 1 on page l4 and the following formula (Ref. 1: pp. 423-424]: 

Оа со 


Table 1. ATMOSPHERIC PARAMETERS AND VALUES 


E ПЕ. 


ди 15614] 0 ЕШ ШЕШЕЙ 2. PSOE | 5 


1.015484E-07 2.21698Е-07 КЕЛЕТ Е-09 


= ae aie 





2. Сопота зе ок риса а во. 
a. m = mass of the satellite. set equal to 100kg. 
b. AR = the cross-sectional area of the vehicle perpendicular to the direction of 
motion. set equal to 20H? 
3. PERTURBING FORCE DUE TO HEAVENLY BODY 
The satellite will experience perturbation forces due to the gravitational effects 
of the sun and the moon. The perturbation force from a perturbing body is the differ- 
ence between the gravitational force due to the perturbing body at the satellite and the 
gravitational force the satellite would expemence if 1t were at the center of the camer 


From Figure 2 on page 13, {Ше решио е E O Шене 


The variable and constants are defined below: 
1. Малаві 
a. r, = distance from the earth center for the perturbing bodv 
b. i, = unit vector from the earth to the perturbing body 


distance from earth center to the satellite 


с. г = 
d. i — unit vector from the earth to the satellite 

2 Constans: 
а. џи, = gravitational constant of the perturbing body = If,G 


The subscript pis to be replaced by Sai the pertureine body is the sun, and by ш Ш 
nerturbing body is the moon, We will assume thay, < <= + then the equation cao 


“р * мс 
С О 5 


Satellite Perlurbing 
Е ни Гу. body 


Earth 





Figure 2. Perturbation forces. 
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The unit vectors i, and р can be written in terms of the ‘IJK’ system as: 
i, = ( cos(£2) cos(u) — sin(Q) cos(z) sin(ug)4 t ( cos(ug) sin(£2) + соз(0)) cos(i) sin(us))J + 


( sin(z) sin(u))K 


< 


i, — (cos(Q.) cos(u,) — sin(Q,) cos(i;) sin(uop)) 7 + 


c 


(соз(иор) sin($2,) + cos(wp) cos(ip) (идр) + ( sin(é,) sin(uop))K 


where £2. 1, and uw are the orbital elements of the satellites and Q, i, and ta, are the 
orbital elements of the perturbing body. The formulas above use the ‘IJK’ system, and 
as such the resultant forces must be transformed to the 'RSW’ system. Models of the 
sun and moon orbits are required to calculate r, and i. The models used in the program 
for the sun and moon’s orbits follows: [Ref. 3: pp. 73-74] 
a. SUN'S POSITION 

In order to model the suns orbit, a number of simplifications had to be 
made in the actual parameters of the suns orbit. First the sun will be assumed to be in 
a circular orbit. This means that the radius (r) to the sun will be constant. and the ec- 


Benuuem (e) will'equabO.0 mstead of its true value of 0.017. The other assumption will 


be to place the sun on the “I axis ofthe IJK system at the beginning of the program 


and have it progress through its orbit as the program runs. These changes will not effect 


the perturbing force in anv noticeable magnitude. 


The following Variables and constants where used in the program to model 


the suns orbit after applying the simplifications: [Ref. 3: pp. 75-78] 


I 


Constants: 


E А (S 
a. Gravitational Constant: G = 6.67£ — 11 ——— 


Ко: 
b. Sun’s Mass: m, = 1.99E30Ag 
c. Sun’s Gravitauonal parameter: 


- 2 
Ут 
! 


x 
ht 
“> 


и; = 1.32733E20 


d Sün s eccentrici 2 0009 

e. Radius of orbit. assume sun is in circular orbit: у, = у | ]т 
f Sun's inclination: si = 23.45 deg. = 4.09279709d-01 radians 
v. Longitude of ascending node: £2. — (0:0 

h. Xreument of perreee; c). = 00 

Ма 0169: 


a. The true anomalv of the sun’s position as a function of the time the satellite has 
been in orbit: 


= о 


Where TT = true time. the time the satellite has been in orbit (sec) 


— 


b. Sun s'Position vectora m — dec sm EIE 


| 


с. Unit Vecror оте ено = 


~ 
~ ° 


b. MOONS POSITION 


In modeling the orbit of the moon. simular assumptions where used as with 


the sun. The moons orbit will be assumed to be circular, actually the eeceentrer wp. 


equal to 0.055. By placing the moon initially on the Г ах of the ‘IJK’ system along 


with the sun, the gravitational forces of the two bodies will combine to a maximum. 


However; since the moons orbital period. 1s only 27.3 davs, the moon will not stay in this 


alignment and the magnitude of the combined forces will vary with ume. The inclination 


of the nroons orbit is not constant. bur drifts between 18.5 and™@s.6 decrees mm tem rea 


Also the longitude of the ascending node (©) oscillates between 13 and -13 degrees. To 
simplify this the inclination will be chosen as a constant 23.5 degrees and the longitude 
of the ascending node as 0.0 degrees. For the time period involved in calculating the 
perturbed orvit, these assumptions will not make any significant difference. 

The following variables and constants were used in the program to model 
the moons orbit, after applving the simplifications: 


PS Onstants: 


ща b (Хот?) 
пи гахнацопа! Сопсзап. б = 6.672 = 11 m 
b. Moon's Mass: m, = 7.55 E22kg 
EET | (Nm?) 
c. Moon's Gravitational Parameter: k, = ОМ, = 4906 12 —— 
kg 


d. Moon's eccentricity: e, = 0.0 
e. Radius of orbit, assume moon is in circular orbit: r, = 3.844 E8km 
f. Moon's inclination: i + 25.5Чег. = 4.10152374E-1 radians 
v. Moon's longitude of ascending node: 6), = 0.0 
h. Moon's argument of perigee: c, = 0.0 
1. Moon's period: T 2 27.3 davs [period] 
Ша сатта р!с<: 


a. The true anomaly of the moon's position as a function of the time the satellite 
э 


mg кш ор Uu imc И 
их Оон) 


-» 


b. Moon's position Vector: y 2 r cos vy,P Fr SIN TinQ 


{ Е ) 
ШИЕ ТІ хесог from earth to moon: 5, = 


pun 


The models of the sun and moons orbit calculates the position vector in the POW’ 





svstem and therefore the position vector must be transformed to the ‘IJK’ system. 


С. RATE-OF-CHANGE OF ORBITAL ELEMENTS 

The derivations and equations of the rates-of-change of the orbital elements are 
contained in reference ] pages 398 to 406. Therefore: only a summary of the actual an- 
alytic expressions for the rate-of-change of the parameters in terms of the perturbations 
Will follow: 


|. Rate-of-chanege of the semi-major aas: 


„ә 


CA 


~ 


/ зе сүү SEI Ee: 
«ча -€ SIN v, Е EN = 
A ALL cepe 3 ыы з. ада ет са arr, Б 
E EXE Е И Aj | 
dl i ` THO 
H NS | ms с 


Where 27’ 1s the mean motion of the satellites orbit. 





Waa 
Rate-of-change Of the eceentniciin: 
T vli e^ sin vg Ве. d eese) 
es [| SF DM — Y МА 
at Ziel мате 


Rate-of-change of the inclination: 


di F COS Lig 


at m poste 2 
Rate-of-change of the longitude of the ascending node: 


ао Y Sin Up Е 
dt = [ ^ ; | | ] и’ 
пы а sin? 


Rate-of-change of the argument of perigee: 


oy (:) 
-- |) --(а---)-- 
A uet met GT. 
MEE 
= | — er COS V 
do N 0 Б 
dt = п ае ІШ 
(E), = (2 sin (1 ЛЕ 
dt | de cos EE 
б —yr CO Í SIN iha 
(r L ЈА, 
n'a V | — е” 


Rate-of-change of the eccentric anomaly: 


: 2 ! de $ ps ` + de (е Жыры \ ' 
Е | [( sin vg 4 = (1 + есоѕ т) — ( соу) + еј y COS o + esin Vo)] 


dt — sin(E.4) moo 


Rate-of-change of the mean anomaly: 


d M Á ЧЕ! de (EA\dEA du | 
IM LG. -c—ndLi-excos———— — -—— (t- 1 
dt al dt dt di 


This equation reduces to the following for circular and ecliptic orbits 


(0< =е2< 1). 

ШЕ po 0-е 1-<: a! 
QC > -- Е € р. EM Y , uH 
С с! |9 и 

(41 H d мае а(1 ЩЕ. е”) dl 


Where the Rate-of-change of the mean motion: 


ам Е E da 


а! m di 





[ref. 1 p. 396-407] 


D. NEW ORBITAL ELEMENTS 

The change of each element is calculated by multiplving the rate-of-change of the 
element bv the time step (DT). The change in the orbital elements are then added to the 
current values of the elements to give the new orbital elements. With the new elements 
calculated, the satellite is stepped forward and the new position and velocity are calcu- 
lated in the same manner as the unperturbed orbit (chapter 3). Also as with the unper- 
turbed orbit, the process is repeated until the satellite is at the perigee point. indicated 


by the time of flight (TF) equal to the period of the perturbed orbit. 
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у. БОСС ВЕ а 


The ability of the student to change the velocity of the satellite at any position in 


the orbit is a vital element in this program. With velocity changes the student can in- 


vestigate the effects of varying the satellites velocity as in transfer orbits and inclination 


changes. In order to simplify the program the unperturbed orbit is used throughout this 


routine. The velocitv change algorithm used in the program follows: 


ЈЕ 


ty 


v», 


I. 


Rotate to velocitv change location. 

The user is given the choice of changing the velocity of the satellite atte 
perigee, apogee or at anv true anomaly. If the user chooses perigee or apogee as 
the change locauons, the true anomaly is set equal to zero or pi radians respect- 
fullv. With the location of the velocity change, the satellite is first stepped around 
to the desired true anomaly. The stepping 15 identical with the unperturbed orbit 
with the exception that the stepping terminates when the true anomaly is greater 
or equal to the desired true anomaly. With a step size of one fiftieth of the period, 
the satellite is actually stepped around to a location near the desired location. This 
variance can be reduced bv decreasing the step size but this would increase the 
computation tme. This error will be a major factor in precise calculations of 
transfer orbits, or anv other orbital maneuver where precise Velocity changes are 
required. However: this program 1s not a tool to calculate precise orbital maneu- 
vers, but rather a learning tool for the student to get a feel for the results of velocity 
changes in a satellites orbit. 


Се Ее 

With the satelite at the desired location, the program calculates and displays 
for the user the satellite's current velocity, escape velocity and circular чејосо 
velocity required to circularize the orbit). The program will not allow velocities 
greater than or equal vo the escape velocity. The user is given the option (Ok -mii 
a new velocity in the IJK’ svstem or to change the magnitude of the velocity inthe 
orbital plane. If the user chooses to change the velocity in the orbital plane me 
program Will prompt the user for the magnitude of the velocity change. and multi- 
plv this change bv a unit vector in the direction of the satellites velocity. DISSE 
locitv change vector 1s then added to the satellites velocity vector. to calculate the 
new velocity vector. 


Calculate new elements. 
The orbital elements are calculated with the new velocity vector and the satel- 
lite’s position vector. 


Complete the orbi 

The program will complete the orbit to the new perigee point using the satel- 
lite s position. new velocity and new elements. There are a number of problems 
that arise if the satellite 1s just stepped around to the perigee point. For example. 
with velocity changes in the orbital plane the apogee and perigee directions can 
physically swap. This is a problem when plotung with the perifocal coordinate 
svstem because the .Y, axis points toward perigee. To avoid problems like this the 
arrays used in plotting the orbit must be cleared and the satellite's current position 


and velocity be treated as initial conditions. Hlowever:; to compare the old and new 
Orbits there is a desire to retain as much of the previous orbit as possible. The 
velocity changes where divided into the following four cases to handle these prob- 
lems: 


4: 


9. 


feater tian the circular velocity sme perigee point will remain the same so the 


Change velocity in the orbital plane at the perigee point with the new velocity 
C 
satellite is stepped around using the unperturbed subroutines. 


“- 
- 


Change velocity in the orbital plane at the perigee point with the new velocity 
less tian or equal to the circular velocits.. The perigee and apogee directions 
pul иер 50 the plotting arrays dre first cleared and stored with the current 
location data. Because the satellite is now at the apogee point the satellite is 
stepped around to the perigee point storing the second half of the orbit. The 
entire next orbit is calculated and stored to get a complete orbit. 


Change velocity in the orbital plane at the apogee point with the new velocity 
less than the circular velocity. The perigee and apogee directions will remain the 
same. so the satellite is stepped around to the perigee point completing the orbit. 


This last case catches all the following velocity changes; velocity change in the 
orbital plane at the apogee point with the new velocity greater than or equal to 
the circular velocity. velocity changes at anv other true anomaly in the orbital 
plane. and anv velocity change out of the orbital plane. The plotting arravs are 
cleared and stored with the current location data. No matter where in the orbit 
the satellite 1s, the satellite is first stepped around to the perigee point. and to 
ensure a complete orbit is plotted the entire next orbit is also calculated and 
БОТЕВ. 


VI. GRAPHICAL PLOTS 


The program provides two types of graphical displays of the orbit, a display in the 
perifocal coordinate svstem and a display of the satellite’s ground track. Each displav 
type is useful in observing different aspects of the orbit. The perifocal display will allow 
the user to see how certain orbital parameters change with different initial positions and 
velocities, and also how the parameters change with velocity changes at varying posi- 
tions in the orbit. The ground track will enable the user to gain an appreciation for the 
physical location of the satellite above the earth, and see how the orbital parameter af- 
fects the path of the satellite. The ground track will also display the precession of a se- 
quence of orbits. Both displavs plot the position steps to give the user an understanding 
of how the satellite speeds up at perigee and slows down around apogee. 

The DISSPLA package on the mainframe computer was used to enable the plotting 
of the orbits. The versatility of plotting subroutines of DISSPLA makes the actual 
programming of the orbit a simple matter of initializing DISSPLA for the type of mon- 
itor being used, setting up the plotting area, initializing the axis and axis scale, and then 
plotting the desired curve from points contained in arrays. This is a simplified explana- 
tion of DISSPLA, but for further details on DISSPLA programming refer 10 ІС 
DISSPL.A user's manual [Ref. 5. DISSPLA also supplies subroutines to draw a variety 
of projections of the world and fill the projections with coast lines, latitude lines and 
longitude lines. There are a couple of DISSPLA requirements that did require special 
handling in the program. The requirement that the data be supplied in arrays forced the 
program to load arrays with the required position and parameters and to keep a counter 
for the number in the arravs. The array format requires the size of the arrav be specified 
in the beginning of the program. The array size needs to be large enough to hold a 
number of orbits, but not so large as to waist storage space. The program will continue 
to add orbital data to the arrays until the user chooses to delete the previous orbits. If 
a new initial position and velocity is entered or if the arrays will overflow with the next 
orbit the arravs will automatically delete all previous orbits. DISSPLA also requires that 
all data be in single precision format. The program calculates all orbits in double pre- 
cision in order to limit the effect of round-off error, but by using the single precision data 


for plotting will not affect the accuracy of the plot in any wavy. 


9” 


The subroutines used to displav the orbits will be covered in the following three 


sections: 


meer ERIFOCAL PLOT 

The plotting of the orbit in the perifocal coordinate svstem is the easier of the two 
tvpes of plots. Since the perifocal coordinate system has the orbital plane as the fun- 
damental plane, the only requirements to describe the orbit in the perifocal coordinate 
system are arrays with the true anomaly and the radius to the satellite. To give the user 
a sense of the size of the plot. the axis length varies with the eccentricity and semi-major 
axis length. Also a plot of the earth 15 plotted to the same scale, with the pole or center 
of the plot on the origin of the axis. The latitude of the earth at the center of the plot 
will vary with the inclination of the orbit. Thuis plot will allow the user to see a relative 
view of the satellite’s coverage in the minus ‘Z’ axis direction of the perifocal coordinate 


system. 


B. GROUND TRACK 

The ground track plot is a very complex subroutine compared with the perifocal 
ни Ресашсе ше сгоџпа track 1s not a continuous curve a procedure to handle the 
satellite ending at one end of the plot and wrapping around to the other end was devel- 
oped. The wrap around problem is avoided in most orbits by plotting the orbit in seg- 
ments with the following two rules. Each segment begins at the beginning of a new plot 
ИАЕА of the plot area. and ending when the satellite would wrap around to the 
other side of the plot. At the beginning of a segment if the position of the satellite ts 
within five degrees of the edge of the plot, that position and anv other positions within 
that five degree boundary will not be plotted. The segment will end when the satellite 
15 within ten degrees of the edge of the plot. The above restricuons imposed on the 
segments of the plot will not substantially affect the interpretation or usefulness of the 
plot. The ground track 1s plotted on top of a cylindrical equidistant projection of the 


world, with the world coast lines and a longitude-latitude grid for reference. 


C. DATA 

Informauon concerning the orbit is displaved on the lower half of the plot. The 
information 1s designed to supply the user with enough of the basic orbital elements and 
other parameters affecung the orbit to be able to evaluate what basic tvpe of orbit the 
satellite 15 in, and the effects of velocity changes and perturbing forces have on the orbit. 


The following data are plotted: inclinauion(1), semi-major axis (a), eccentricity (¢), period 
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pengee velocity and radius, average time rate-of-change of orbital 
с f 


agniiude of perturbing forces per unit mass. 


VH. CONCLUSIONS AND RECOMMENDATIONS 


The program supplies the student with an interactive tool to study the orbital mo- 
tion of satellites around the earth. The student can investigate a variety of orbits by 
Varving the orbital parameters, command velocity changes, and observe the effects of 
perturbing forces. 

The student is provided with two options for entering the initial position and veloc- 
itv of the satellite. The program could be expanded to provide the student with the ad- 
ditional options of entering either orbital parameters or a ground observation data and 
have the program calculate the initial position and velocity from this data. Also the 
student 1s limited to orbits with eccentricities less than one (elliptic orbits). The program 
could be also be expanded to include more eccentric orbit for Lunar, interplanetary. and 
missile trajectories. The perturbing orbit is calculated for orbits around the earth with 
relatively small perturbing forces in relation to the earths gravitational force. This fact 
will cause the program to produce false results if the student tries to calculate lunar 
Meajectories. Special routines would have to be emploved when the perturbing force (the 
moons gravitational attraction) 1s comparable to the earths gravitational attraction. 
Mise will not become a factor for studying current satellite orbits out to th 
geosvnchronous radius of 42241.1km. 

Tne velocity change subroutines move the satellite to a location close to the desired 
(Поп before a velocity change is imposed. Bv reducing the step size in the velocity 
change subroutine. this error could be reduced. Precise orbital transfer maneuvers can 
be modeled by reducing this error caused by the positioning of the satellite prior to 
changing the velocity. The program will currently provide the student with useful plots 
for gaining eXperience with various transfer orbits bv varving the magnitude and location 
of the velocity changes. 

The output of the calculations of the orbit are arravs loaded with the satellite’s po- 
sition and select orbital parameters. The DISSPLA subroutines that plot the points are 
not unique. The program would become portable to personal computers with these 
graphics subroutines written in FORTRAN and included in the program. 

A final recommendation is that the display of the ground track could be modified 
to show ground coverage, number of satellites in a constellation, and other elements 


necessary for planning a real-world artificial satellite application. 


APPENDIX A. ORBIT PROGRAM 


PROGRAM ORBIT 
THIS PROGRAM IS AN INTERACTIVE TIME STEP SIMULATION OF 


SATELLITES AROUND THE EARTH. 
ARE CALCULATED AND PLOTTED. 


PERTURBED AND UNPERTURBED ORBITS 
VELOCITY CHANGES ARE ALSO PERMITTED 


AT SPECIFIED TRUE ANOMALIES. 


> 
с 
IA 
on 
= 


пп ни ни пи и и ни и ни пио ни 


H Hu н пи ни п ни нп и ни пп ни ни ни 


= OPTIONS: 


OF VARIABLES USED BY THE MAIN PROGRAM FOLLOWS: 
SEMI-MAJOR AXIS 

ARGUMENT OF LONGITUDE 

ARGUMENT OF PERIGEE 

VELOCITY CHANGE LOCATION TRUE ANOMALY 
TIE Sore 

ECCENTRICITY 

ECCENTRIC ANOMALY 

I VECTOR OF ECCENTRICITY 

J VECTOR OF ECCENTRICITY 

K VECTOR OF ECCENTRICITY 

R VECTOR OF TOTAL FORCE 

S VECTOR OF TOTAL FORCE 

W VECTOR OF TOTAL FORCE 

ANGULAR MOMENTUM 

I VECTOR OF ANGULAR MOMENTUM 

J VECTOR OF ANGULAR MOMENTUM 

K VECTOR OF ANGULAR MOMENTUM 
INCLINATION 

PERTURBED OR UNPERTURBED OPTION 
PLOT NEXT ORBIT, CHANGE INITIAL VALUES, 
VELOCITY, PLOT ANOTHER VIEW OF ORBIT, QUIT 
LONGITUDE OF ASCENDING NODE 
LONGITUDE OF PERIGEE 

MEAN ANOMALY 

MEAN MOTION 

GRAVITATIONAL PARAMETER 

ASCENDING NODE 

I VECTOR OF ASCENDING NODE 

J VECTOR OF ASCENDING NODE 

K VECTOR OF ASCENDING NODE 

STEP COUNTER 

SEMI-LATUS RECTUM 

РЕКТӘВ OF СЕБЕ 

ЕІ 

RADIUS OF APOGEE 

RADIUS OF EARTH 

ORBITAL RADIUS 

I VECTOR OF ORBITAL RADIUS 

J VECTOR OF ORBITAL RADIUS 

K VECTOR OF ORBITAL RADIUS 

TIME COUNTER IN ORBIT 

TREE ANOMALY 

TOTAL CHANGE IN SEMI-MAJOR AXIS 
TOTAL CHANGE CINCARGOUIEST OPDOPERTGEE 


ОКВ00010 
ORB00020 

ОКВ00030 
ОКВО0040 
ОКВО0050 
ОКВО0060 
ОКВ00070 
ОКВ00080 
ОКВО0090 
ОКВ00100 
ОКВ00110 
ОКВО0120 
ОКВ00130 
ORB00140 
ORB00150 
ORB00160 
ОКВ00170 
ORB00180 
ОКВ00190 
ОКВ00200 
ОКВО0210 
0КВ00220 
ОКВО0230 
ORB00240 
ORB00250 
ОКВО0260 
ОКВО0270 
ОКВО0280 
ОКВО0290 
ОКВ00300 
ОКВ00310 
ОКВ00320 
ОКВ00330 
ОКВ00340 
ОКВ00350 
ОКВ00360 
ОКВ00370 
ОКВ00380 
ОКВ00390 
ОКВО0400 
ОКВОО04 10 
ОКВО0420 
ОКВОО04 30 
ORB00440 
ОКВ00450 
ОКВО0460 
ОКВОО04 70 
ОКВ00480 
ОКВО0490 
ОКВ00500 
ОКВ00510 


c 


[DE = TOTAL CHANGE, IN ECCENTRICITY 

TDH = TOTAL CHANGE IN ANGULAR MOMENTUM 

TDI = TOTAL CHANGE IN INCLINATION 

TOMA = TOTAL CHANGE IN MEAN ANOMALY 

В = Tere CHANGE. IN MEAN HOLION 

TDLAN= TOTAL CHANGE IN LONGITUDE OF ASCENDING NODE 
ШЕ = ТЕ ОЕЗЕРТОНТ 

TrDRA= TOTAL FORCE OF DRAG 

TFEA = TOTAL FORCE OF EARTH'S OBLATENESS 

TFMO = TOTAL FORCE FROM MOON 

ВЕ-О = TOTAL. FORCE FROM SUN 

TL = TRUE Longitude AT EPOCH 

TI = TRUE TINE SINCE SATELLITE HAS BEEN IN ORBIT 
V = SATELLITE VELOCITY 

VI omm elk OF SATELLITE VELOCITY 

VI = J VECTOR OF SATELLITE VELOCITY 

УК = К VECTOR OF SATELLITE VELOCITY 


A LIST OF THE ARRAYS USED FOLLOWS: 


AINRAY = INCLINATION 

APRAY = ARGUMENT OF PERIGEE 

RARAY = RADIUS 

RIRAY = I VECTOR OF RADIUS 

RJRAY = J VECTOR OF RADIUS 

RKRAY = K VECTOR OF RADIUS 

TARAY = TRUE ANOMALY 

TIMRAY = TIME 

A LIST OF SUBROUTINES CALLED BY THE MAIN PROGRAM WILL FOLLOW: 
CALCEL = CALCULATES THE ORBITAL ELEMENTS 

ШІСУЕШ = ALLOW THE USER T@ CHANGE THE VELOCITY OF THE SATELLITE 
INPUTS = PROMPTS USER FOR INITIAL POSITION AND VELOCITY 

ШЕ С = INITIALIZES THE SUMS IN THE ARRAYS 

NEWELT = CALCULATE NEW ORBITAL ELEMENTS FROM TIME STEP 
NEWPOS = CALCULATE NEW POSITION VECTOR 

NEWVEL = CALCULATE NEW VELOCITY VECTOR 

OPTION = GIVE THE USER THE OPTIONS Permitted IN THE PROGRAM 
HIS = PLOTS THE OKBITS 

PRETUR = CALCULATES THE PERTURBED ORBIT 

STORE = STORE THE POSITION DATA IN ARRAYS 

UNPRET = CALCULATE THE UNPERTURBED ORBIT 


BEGIN MAIN PROGRAM 

DOUBLE PRECISION PI,MU,RI,RJ,RK,R,VI,VJ,VK,V, HI, HJ, HK,H, 
+ NIQNJQNE,N,PB,EI,EJ,ER,E,A, I, LAN, AP, TA, AL, LP, TL, PER,EA, 
+ MM,MA,T,DT,TF, FR,FS,FW,TT,CHTA,RA,VAÀ, TEMPTA, RE 


DIMENSION TARAY(500) , RARAY(500) , RIRAY(500) , RJIRAY(500) ,RKRAY(500), 
+ AINRAY(500) ,APRAY(500) , TIMRAY(500) 


CHARACTER*1, LOOP, YORN , ORLOOP 


РІ -735/1515925395099035 


08800520 
05800530 
08800540 
ОКВО0550 
ОКВО0560 
ОКВО0570 
ОКВО0580 
ОКВ00590 
ОКВ00600 
ОКВ00610 
ОКВ00620 
ОКВ00630 
ОКВ00640 
ОКВ00650 
ОКВ00660 
ОКВО0670 
ОКВО0680 
ОКВ00690 
ОКВ00700 
ОКВ00710 
ОКВ00720 
ОКВ00730 
ОКВ00740 
ОКВ00750 
ОКВ00760 
ОКВ00770 
ОКВ00780 
ОКВ00790 
ОКВ00800 
ОКВ00810 
ORB00820 
ОКВ00830 
ORB00840 
ОКВ00850 
ОКВ00860 
ОКВ00870 
ОКВ00880 
0ОКВ00890 
ОКВ00900 
ОКВ00910 
ОКВ00920 
ОКВ00930 
ОКВО0940 
ОКВО0950 
0ОКВ00960 
ОКВ00970 
ОКВ00980 
0КВ00990 
ОКВ01000 
ОКВ01010 
ОКВ01020 
ОКВ01030 
ОКВО1040 
ORBO1050 
ОКВ01060 
ОКВ01070 


10 


20 


=o, 7660120705 
= 6. 378145D+03 


INTRO TO PROGRAM 
INTRO 


ENTERED MAIN PROGRAM LOOP 
LOOP = 'Y' 


ТЕ 


(LOOP .EQ. 'Y') THEN 

INITIALIZE STEP COUNTER AND TRUE TIME 
NUM = 1 

TT = 0.0 


PROMPT USER FOR INITIAL POSITION AND VELOCITY 
CALL INPUTS(RI,RJ,RK,R,VI,VJ, VK, V, MU, LOOP, PI) 


EXIT PROGRAM 
IF (LOOP THEN: 

GOTO 10 
ENDIF 


'N!) THEN 


CALCULATE AND STORE ORBITAL ELEMENTS 
CALL CALCELCRI, RJ, RK,;R ТТ, УТ, УК, ЕВЕ ВЕ ВТ. ВАМ, 
LP, TA, PER, EA, MA, AP, AD IFP FLNO MEAN H, RTE 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY, TARAY, 
NUM,I,AP,AINRAY,APRAY,TT,TIMRAY) 


PRINT DATESEOR 05БЕ ТО ВЕУТЕН 
ER 


PRINTS, ИЕ ИВ" 
PRINT", VJ = WI NDS 
PRINT, VK = , ҮК : ШІСІ 
PRINT*,' V s', V,' KM/S' 
PRINT*,'RI z', RI,' KM' 
PRINT*, RJ = . ЕТЕ 
PRINT, RK ШЕК, ee 
PRINT,” R=, Rk,’ KM 


PRINT*, ‘ECCENTRICITY =' ,E 

DEGI = SNGL( (180. 0/PI)*I) 

PRINT, ‘INCLINATION =' ,DEGI,' DEGREES' 
PERHRS 2 SNGL(PER/3600. 0) 
PRINT*,'PERIOD z',PERHRS,' HOURS' 
PRINT*,'ARE THESE VALUES CORRECT? 
READ* , YORN 

CALL, EXCMS@ CERSCEN ) 


EMR Yo Ола = 


IF (.NOT. YORN .EQ. 'Y') THEN 
GOTO 20 
ENDIF 
CALCULATE TIME STEP AND SET TIMER TO ONE TIME STEP 
DT = PER/50 
Tm 


ӘТЕР БАТЕБ ШЕ ІСЕРЕЕКТБЕР РОСА TASOFRECORD 
Ir ACTA GT. OOE STA D АІ Б СІНЕ 
EI ЖІ 


ОКВО 1080 
ОКВО1090 
ORBO1100 
ОКВ01110 
ОКВ01120 
ОКВ01130 
ORBO1140 
ОКВ01150 
ORBO1160 
ORBO1170 
ОКВО1180 
ОКВ01190 
ОКВО1200 
ОКВО1210 
0КВО1220 
ОКВ01230 
ОКВ01240 
ORBO1250 
0КВ01260 
ОКВ01270 
0КВО 1280 
0КВ01290 
ОКВ01300 
ОКВ01310 
ОКВ01320 
ОКВ01330 
ОКВО1340 
ОКВ01350 
ОКВО1360 
ORBO1370 
ОКВО1380 
ОКВО1390 
08801400 
ОКВО1410 
ORBO1420 
ОКВО 1430 
ОКВО 14540 
0КВО1450 
0КВО1460 
ORB01470 
ОКВО 1480 
ОКВО 1490 
ОКВ01500 
ОКВ01510 
ОКВ01520 
ОКВ01530 
ОКВ01540 
ORBO1550 
ОКВ01560 
ОКВ01570 
ОКВ01580 
ОКВ01590 
ОКВ01600 
ОКВ01610 
ОКВ01620 
ORB01630 


зе 


СОО БИ БЕО в А, ТА, ТЕ, ЭТ, РГ. РЕВ) 

ООШ РОК, В. ВЕ. в АХ. АР,Т, ТА, А.В) 

И МЕ, РВ РУ, УК) 

ем 1 

САПЫ SIORE( RK] ,RI,RR GR, TA, RIRAY , RJRAY , RKRAY , RARAY , TARAY, 
mum, 1, AP ,APNRAY ,APRAY, TT ,TIMRAY ) 

ПЕРЕ ПТ 

GOTO 350 


ENDIF 


CALCULATES ELEMENTS FROM PERIGEE POINT 
Cee в АВВ. ВА. ВУ, VK V EL EJ EK, Eâ, I; LAN, 


ESTA FERK ESNIA AP AL ГОРУР, МО МММ, ННІ, Н) 


DT = РЕК/50 
T = DT 


STORE FIRST Unperturbed ORBIT 
ВИ PRET DI PER AL. LAN, AP. I КЕ КЈЈЕК К,УТ, УЈ,УК,У, 
в Тр. 11 А. ВА ТЕ. NUM, RIRAY,RJRAY, 


RKRAY , RARAY , TARAY , AINRAY , APRAY , TIMRAY , TT) 


INITIALIZE SUMS FOR FORCE AND ORBITAL ELEMENT CHANGES TO ZERO 
(РОЛАН А ИУ РОВ АИ У О ОРАО TFDRA, IDI, TDA, TDE, TDM, TDMA, TDLAN, 


ТОН , TDAP ) 


ГЕШ ЕТКЕ UNPERTURBED ORBIT 
го при ПОТ (RIRA RIRAY RERA КАВАУ ТАКА ХОМО РТ, 1 ТР, А,„Е,„ТЕ, 


ОКВО1640 
ОКВО1650 
0К801660 
ОКВО 1670 
ORBO1680 
ОКВО 1690 
ОКВО1700 
ORBO1710 
ORBO1720 
ОКВО1730 
ОКВО1740 
ОКВ01750 
ОКВ01760 
ОКВ01770 
ОКВ01780 
ORB01790 
ОКВО1800 
ОКВ01810 
ОКВО1820 
ОКВО1830 
ORB01840 
ОКВ01850 
ОКВ01860 
ОКВ01870 
ОКВ01880 
ОКВ01890 
ОКВО1900 
ОКВ01910 


ATNESPOCADRATOIRHIRAYOSPBSA.DESDSTEMO,TFDRA,PER,IDI;TDA;ORBOI920 


DUE UBL TDAASTDLAN, TDH, TDAP,MM, MA, LAN, H, AP,R, V) 


CHANGE INITIAL VALUES 
CHANGE VELOCITY AT A SPECIFIC TRUE Anomaly 
PLOT ANOTHER VIEW OF SAME ORBIT 


BEGIN NEW ORBIT OPTIONS 
еи Unperturbed ORBIT 
= Perturbed ORBIT 
= QUIT 
IOPT2 = PLOT Nee ORBIT 


I WN KR Оз 


ALSO ASKED IF WANT TO CLEAR ALL PREVIOUS ORBITS 


CALCULATE ELEMENTS AT PERIGEE 
CALL CALCEL(RI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E, A, I, LAN, 


LP,TA,PER,EA,MA, AP, AL, TF, P, PI, MU, MM, N, H, HI, HJ) 


CHECK FOR POSSIBLE ARRAY OVERFLOW 


IF (NUM .GT. 425) THEN 
PRINT*,'ARRAYS ARE FULL' 
PRINT*,'PREVIOUS ORBITS WILL BE ERASED!' 
NUM 1 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY, RKRAY , RARAY,TARAY, 
NUM,1,AP,AINRAY, APRAY, TT, TIMRAY) 


ESD EE 


БЕРЕКЕ КЕНЕ ТЕЕВОРЛТОТОА 


ОКВ01930 
ОКВО1940 
ОКВ01950 
ОКВО1960 
ORBO1970 
ОКВО1980 
ORB01990 
ОКВО2000 
ОКВ02010 
ORBO2020 
ОКВО2030 
ОКВ02040 
ORB02050 
ОКВО2060 
ОКВ02070 
ОКВО2080 
ОКВО2090 
ОКВ02100 
ОКВ02110 
ОКВ02120 
ОКВ02130 
ORBO2140 
ОКВ02150 
ОКВ02160 
ОКВ02170 
ОКВО2180 
ORBO2190 


зе 


ze 


+ + + 


+ + + + 


+ + + + + 


+ + + + 


+ + + + + + 


+ 


CALL OPTIONCIOPT1, IOPT2 Ом ЕТЕД И ККЕ ТК 
TARAY , AINRAY , APRAY , TIMRAY ) 


Initialize SUMS FOR FORCE AND ORBITAL ELEMEN! CHANGES TOT EER 
CALL INTSUM(TFEA,TESU,TEFMO, TFDRA,TDI,TDASTDE S TDMM, TDJA ; TOP И 
TDH, TDAP) 


SET TIME COUNTER TIO ОЕ НЕ IET 
шъ 


OPIION: PLOT THE УВА ТОКОВИ 
ТЕ (ТОРТ? . EQ. ЈЕНЕ 


CALCULATE AND PLOT UNPERTURBED ORBIT 
IF(IOPT1 .EQ. 1) THEN 

CALL UNPRET(DT,PER,AL,LAN,AP, I, RI,RJ, 
RK,R,VI,VJ,VK,V,MU,PI,H,A, 

E,N,TA,P, MM,MA,EA, TF, T, NUM, RIRAY,RJRAY, RKRAY, 
RARAY , TARAY , AINRAY , APRAY , TIMRAY , TT) 

CALL PLOTS(RIRAY,RJRAY,RKRAY , RARAY , TARAY , NUM, 
PI,I,LP,A,E,TF,AINRAY , APRAY, TIMRAY, 
TFEA,TFSU,TFMO,TFDRA,PER, 

TDI, TDA, TDE,TDMM, TDMA, TDLAN, TDH, TDAP, 
MM,MA,LAN,H,AP,R,V) 


CALCULATE AND PLOT PERTURBED ORBIT 
ELSEIF(IOPT1 .EQ. 2) THEN 
CALL PRETUR(DT,PER,AL,LAN,AP,I, 
RI,RJ,RK,R,VI,VJ,VK,V,FR,FS,FW, 
MU,PI,H,A,E,N, TA, P, MM, MA,EA, TF, T, NUM, 
RIRAY,RJRAY,RKRAY , RARAY , TARAY , AINRAY , APRAY, 
TIMRAY,TT,TFEA,TFSU,TFMO, TFDRA, 
TDI,TDA,TDE, TDMM, TDMA, TDLAN, TDH, TDAP) 
CALL PLOTS(RIRAY,RJRAY,RKRAY,RARAY, TARAY , NUM, 
PI,I,LP,A,E,TF,AINRAY, APRAY, TIMRAY, 
TFEA,TFSU,TFMO,TFDRA, PER, 
TDI ,TDA,TDE,TDMM, TDMA, TDLAN, TDH, TDAP, 
MM,MA,LAN,H,AP,R,V) 
ENDIF 


GOTO THE BEGINNING OF THE PROGRAM TO CHANGE THE INITIAL VALUES 
ELSEIF (TOPT2 -EOS ETEN 
ШЕЛ) 


CHANGE VELOCITY AT A SPECIFIC TRUE ANOMALY AND 
PLOT THE NEW ORBIT 
ELSEIF (IOPT2 .EQ. 3) THEN 
CALL CHGVEL(DT,PER,AL, LAN, AP,I,RI,RJ,RK,R, 
VI,VJ,VK,V,MU,PI, 
H,A,E,N, TA, P,MM,MA,EA, TF, T, NUM, RIRAY, 
RJRAY , RKRAY , RARAY , TARAY , AINRAY , APRAY, 
TIMRAY , TT,EI,EJ,EK,LP,HI,HJ,IOPT1, 
TFEA,TFSU,TFMO,TFDRA,TDI,TDA,TDE,TDMM, 
TDMA, TDLAN, TDH , TDAP) 
CALL PLOTS(RIRAY , RJIRAY , RKRAY , RARAY , TARAY, NUM, 
PI ,I,LP,A,E,TF,AINRAY , APRAY , TIMRAY, 


30 


0КВО2200 
ORB02210 
ORBO2220 
ORBO2230 
ORBO2240 
0КВ02250 
0КВ02260 
ОКВО2270 
ORB02280 
0КВ02290 
ОКВ02300 
0КВ02310 
0КВ02320 
ОКВО2330 
ОКВО2340 
ORB02350 
ОКВО2360 
ORB02370 
0КВ02380 
0КВО2390 
ОКВО2400 
0КВ02410 
ORB02420 
ORBO2430 
ORBO2440 
0КВ02450 
0КВ02460 
ОКВО2470 
0КВ02480 
ORBO2490 
ОКВО2500 
ОКВО2510 
ORBO2520 
ORBO2530 
ОКВО2540 
ОКВО 2550 
ORBO2560 
ORBO2570 
ОКВО2580 
ОКВО2590 
ORBO2600 
ORB02610 
0КВО2620 
0КВО2630 
ОКВО2640 
ОКВО 2650 
ОКВО 2660 
0КВО2670 
0КВ02680 
0КВ02690 
ОКВО2700 
ORB02710 
ORB02720 
ОКВО2730 
СКВО2740 
ОКВО2 750 


ve 


90 


+ + 


+ + + + 


ТРЕК ТЕН. ГЕБКА, РЕК, 
ВВ в, Ты, 165. ТОБА ПОН, ТАР, 
ММ, МА, САМ, Н,АР,К,У) 


LOT ANOTHER VIEW OF THE SAME ORBIT 
ESSEN Tren "a THEN 
CALL PLOTS(RIRAY,RJRAY,RKRAY,RARAY , TARAY , NUM, 
РА ИРА ИВ АТИКА ДРКАХ ЛИКА, 
ПРЕД ШрЪре ЛЕНО ТОРА, РЕВ, 
О РОА О ее ТА ПЕТА ТОН ТВАР, 
ММ,МА,ТАМ,Н,АР) 


STOP THE PROGRAM 

ELSEIF (IOPT2 .EQ. 5) THEN 
GOTO 90 

ELSE 
PRINT*, ‘INVALID ENTRY! ' 
GOTO 80 

ENDIF 


CHECK IF SATELLITE Impacted THE EARTH AND GO TO THE BEGINNING 
IF (R .LE. 6450.0) THEN 
PRINT*,'SATELLITE WILL IMPACT THE EARTH!!! ' 
PRINT*,'PROGRAM WILL RESET TO THE BEGINNING' 
GOTO 20 
ENDIF 


НЕСТОР ОЕ ТНЕ ОРТЕОМ LOOP 
GOTO 80 


GIVE THE USER A CHANCE TO RECOVER THE PROGRAM 
PRINT*,' THIS IS YOUR LAST CHANCE! ' 
PRINT*,'DO YOU WANT TO CONTINUE? ' 
PRINT*,' AND GOTO THE Beginning OF THE PROGRAM? ' 
Ека то ENGERCO Y? OR TN s’ 
READ* , LOOP 
PRINT™, LOOP 
GOTO 10 
ENDIF 


DISSPLA SUBROUTINE TO TELL GRAPHICS TERMINAL PLOTTING 
SESolON IS DONE 

CALL DONEPL 

STOP 

END 


ХСХ ХК У ev eye ye 


SUBROUTINE INTRO 
THIS SUBROUTINE WILL GIVE THE USER A Brief INTRODUCTION OF THE 
USES OF THE PROGRAM 


PRINT*,'THIS PROGRAM IS A GRAPHICS DISPLAY OF Satellite ORBITS. ' 
PRINT*,'YOU WILL BE ASKED TO INPUT THE INITIAL VELOCITY AND' 
PRINT*, POSITION VECTORS OF THE Satellite. THE PROGRAM WILL ' 
PRINT*, ‘THEN CALCULATE THE ORBITAL PARAMETERS AND THE ' 


| 


ОКВО2760 
ОКВ02770 
ОКВО2780 
0ОКВО2790 
ОКВО2800 
ОКВО 2810 
ОКВО2820 
ОКВО 2830 
ОКВО2840 
ОКВО2850 
ОКВО2860 
ОКВО2870 
ОКВО2880 
ОКВО2890 
ОКВО 2900 
ORBO2910 
ORBO2920 
ОКВО2930 
ОКВО2940 
ORB02950 
ORBO2960 
ОКВО29 70 
ОКВО2980 
ОКВО2990 
ОКВ03000 
ОКВ03010 
ОКВ03020 
ОКВ03030 
ОКВ03040 
ОКВ03050 
ОКВ03060 
ОКВ03070 
ОКВ03080 
ОКВ03090 
ОКВ03100 
ОКВ03110 
ОКВ03120 
ОКВ03130 
ОКВ03140 
ОКВ03150 
ОКВ03160 
ОКВ03170 
ОКВ03180 
ОКВ03190 
ОКВ03200 
ОКВ03210 
ОКВ03220 
ОКВ03230 
ОКВ03240 
ОКВ03250 
ОКВ03260 
ОКВ03270 
ОКВ03280 
ORB03290 
ОКВ03300 
ОКВ03310 


PRINT*,'Unperturbed ORBIT. 


THE USER WILL THEN HAVE THE' 


PRINT, "CHOIGE OF DISERe ds 


AEN. а e о РТР 1 
м» | N | ` 
2... 1 T 


PRINT*,'THE 
PRINT*,'GRAPHICAL OUTPUT. ' 


PRINT*,' 


1 


-PERIFOCAL (SHOWS RELATIVE SIZE OF ORBIT)' 

-Equatorial (SHOWS ORBIT INCLINED, USER INPUT’ 
LONGITUDE TO VIEW AT)' 

-GROUND TRACK' 


USER IS THEN ASKED TO CHOOSE ONE OF THE FOLLOWING: ' 
-Unperturbed ORBITS' 

-Perturbed ORBITS' 

-VELOCITY CHANGES ' 

USER'S CHOICE WILL BE USED IN DEVELOPING THE! 


PRINT*,'THE USER IS THEN GIVEN THE FOLLOWING CHOICES: ' 


PRINT," 
PRINT*, ' 
PRINT*,' 
PRINT*,' 
RETURN 
END 


-CLEAR ALL THE PREVIOUS ORBITS' 

-CHANGE THE INITIAL PARAMETERS' 

-CHANGE VELOCITY AT A SPECIFIC TRUE Anomaly' 
-PLOT ANOTHER VIEW OF THE SAME ORBIT' 


Wes esee ese edese de desee eese dee de dede ee yede dye edes e edes des edede dede dese eee sede dede ese dede dede dese dese Jededes 


03 


105 


ЩЕ 


+ 


SUBROUTINE OPTION( IOPT1,1TOPT2,NUM,RIRAY , RIRAY , RKRAY , RARAY , 


THIS SUBROUTINE GIVES THE USER A CHOICE OF OPERATIONS THAT CAN BE 


TARAY , AINRAY , APRAY , TIMRAY ) 


PERFORMED ON THE PROGRAM AND RETURNS THE USERS CHOICE WITH 
VARIABLES IOPT1 AND IOPT2 


DIMENSION RIRAY(500) ,RJRAY( 500) ,RKRAY(500) , RARAY( 500) , TARAY( 500), 


AINRAY( 500) , APRAY( 500) , TIMRAY( 500) 
CHARACTER*1,YORN 


ТӘРТІ = 


PROMPT USER FORK OF PIO; 
PRINT*,'WHICH OF THE FOLLOWING OPTIONS WOULD YOU LIKE: ' 


PRINT: ЫШ 
PRINT 
PRIE 
PRINT*, ' 
PRII 
PRINT 
PRINT*,' 


l 


2 


З. 


4 


2 


. CALCULATE THE NEXT ORBIT USING THE SAME' 
PARAMETERS' 

. -CHANGE THE INITIAL PARAMETERS OF THE ORBIT' 

-CHANGE THE VELOCITY AT A POINT IN THE ORBIT' 
(THE UNPERTURBED ORBIT WILL BE USED)' 

. -PLOT ANOTHER VIEW OF THE ORBIT(S)' 

-QUIT' 


РАМТ, БВ. о 

READ* , IOPT2 

PRINT , IOPT2 

CALL EXCMS( ' CLRSCRN'! ) 

IF ( IOPT2 . CGT. 5) THEN 
GOTO 103 


ENDIF 


Prompt USER РОК ТУРЕ ОБ ОКВ ШЕ ЗЛЕ 


ТЕ СЕ: 


12) THEN 


PRINT*,'WHICH TYPE OF ORBIT WOULD YOU LIKE TO SEE, ' 
PRINT*,' 


1. -Unperturbed ORBITS' 


ORB03320 
ОКВОЗЗЗО 
ORB03340 
ОКВ03350 
ОКВ03360 
ОКВ03370 
ОКВ03380 
ОКВ03390 
ORB03400 
ORB03410 
ORB03420 
ORB03430 
ORB03440 
0КВ03450 
ОКВ03460 
0КВ03470 
ORB03480 
ORB03490 
ORB03500 
ORB035 10 
ORB03520 
ОКВ03530 
ORB03540 
ORB03550 
ОКВОЗ5 60 
ORB03570 
ORB03580 
ОКВ03590 
ОКВ03600 
ОКВ03610 
ORB03620 
ОКВОЗ6ЗО 
ОКВ03640 
ОКВ03650 
ОКВ03660 
ОКВ03670 
ОКВ03680 
ОКВ03690 
ОКВ03700 
ОКВ03710 
ОКВ03720 
ОКВ03730 
ОКВ03740 
ОКВ03750 
ОКВ03760 
ОКВ03770 
ОКВ03780 
ОКВ03790 
ОКВ03800 
ОКВ03810 
ORB03820 
ORB03830 
ORB03840 
ОКВ03850 
ORB03860 
ORB03870 


PRINT*, ' 2. -Perturbed ORBITS' 

ВЕР ТЕ МЕБ 1 OR 2:' 

READ TOPTI 

PRINT*,IOPT1 

Сари COO CLRSCRN ) 

ГЕ TOPTI E Те. дк, СТОРТИ КЕ. 

PRINT*,' INVALID ENTRY! ' 
GOTO 105 

ENDIF 

ENDIF 


2)) THEN 


* PROMPT USER TO CLEAR PREVIOUS ORBITS 
107 IF ((IOPT2 .EQ. 1) .OR. (IOPT2 . EQ. 3)) THEN 
PRINT*,'DO YOU WANT TO CLEAR THE PREVIOUS ORBITS?' 
PRINT*,'ENTER "Y" OR "N' : 
READ* , YORN 
PRINT“, YORN 
(ND EXONSC Olrscrn ) 


IF (YORN .EQ. 'Y') THEN 
RIRAY(1) = RIRAY( NUM) 
RJRAY(1) = RJRAY(NUM) 
RKRAY(1) = RKRAY(NUM) 
RARAY(1) = RARAY(NUM) 
ТАКАҮ( 1) S TARAYCNUN) 
AINRAY(1) = AINRAY(NUM) 
APRAY(1) = APRAY(NUM) 
TIMRAY( 1) = TIMRAY( NUM) 
1 

ELSEIF (YORN .NE. 'N') THEN 
PRESE INVALID ENTRY! ! * 
PRINT*,'ALL INPUTS MUST BE CAPITOL LETTERS' 
GOTO 107 

ENDIF 


ENDIF 


* CHECK FOR INVALID OPTION 
ШЕН СТОРТ2 ХЕ. 1). AND. 
+ ШЕН) МЕ а). Акре (порта 

PRINT*, ‘INVALID ENTRY! ' 
GOTO 103 

ENDIF 

RETURN 

END 


CUO Re Eom SAND. 
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Je sev'ezeses e esee eve dese eese ee ese ese dee ees eese ses eese ede dese esee dedeyeseyes ee eee ese ee eve ye e eee 


* COORDINAIE RA EU MS 
warrii У УУ У Је КУ У У Те У је ет 


SUBROUTINE PQWIJKCLAN,AP,INC,P,Q,W,1I,J,K) 
% THIS SUBROUTINE TRANSFORMS PQW COORDINATES TO IJK COORDINATES 


ШЫ КЕШ TON INC SP,O,W,1,J,K,R1I1,R12,R13,R21,R22 ,R23, 
+ В ВЕ, А, АР 


R11 - DCOS(LAN)*DCOS(AP) - DSINCLAN)*DSIN(AP)*DCOS(INC) 
Do РОНА -ОБЈАСАР) = DSINCLAN)*DCOSCAPO*DCOS(INC) 
ЕЕ bom. els.) DSINCINC) 


33 


CTO RI a. NE 3 Y ca xD: 


ORBO3880 
ORB03890 
ORBO3900 
ОКВ03910 
ORBO3920 
ОКВ03930 
ОКВ03940 
ОКВОЗ950 
ORBO3960 
ОКВОЗ9 70 
ORB03980 
ORBOS290 
ОКВО4000 
ORBO4010 
ORBO4020 
ORB04030 
ОКВ04040 
ОКВО4050 
ОКВ04060 
ОКВО4070 
ОКВ04080 
ОКВО4090 
ОКВО4100 
ORBO4110 
ORB04120 
0Е804130 
ОКВО4 140 
ОКВО4 150 
ОКВО4 160 
ОКВО4 170 
ОКВО4 180 
ОКВО4 190 
ОКВ04200 
ORB04210 
0КВО4220 
ОКВ04230 
ОКВ04240 
ОКВ04250 
ОКВО4 260 
ОКВО4270 
ОКВО4 280 
ОКВО4 290 
ORB04300 
ORB04310 
ОКВО4320 
ОКВ04330 
ORB04340 
ORB04350 
ORB04360 
0КВ08370 
ORBO043 80 
ОКВ04390 
ORBO4400 
ОКВ04410 
ORBO4420 
ORBOA4430 


C Yee ee de ve eve Fe Fe Fe He Fede Ns ee Tove Fe Ne Ve eve e Tee Fe He He He Fe Fe Fete He He Tee Fe Ke He He TET AN Te He KA HK HK eK CK 


CO Co Co [2 F2 to 


E CIE Lx d x c 
Проф rD9 Lr 


RET 
END 


SUB 


THIS SUBROUTINE TRANSFORMS IJK COORDINATES TO PQW COORDINATES 


DOUBLE РКЕСІЗІОХ. ТХСУОГЗ,К,ІРИРМИК Е 7 
ИО ЕЗ2, КЭЗГЕ 

R11 » DCOS(LAN)*DCOS(AP) - DSIN(LAN)*DSIN(AP)*DCOS( INC) 
R21 - -DCOSCLAN)*DSIN(AP) - DSIN(CLAN)*DCOS(AP)*DCOS( INC) 
R31 = DSIN(CLAN)*DSINC INC) 

R12 = DSINCLAN)*DCOS(AP) + DCOS(LAN)*DSIN(AP)*DCOS( INC) 
R22 = -DSIN(LAN)*DSIN(AP) + DCOS(LAN)*DCOS(AP)*DCOS(INC) 
R32 = -DCOS(LAN)*DSINC INC) 

R13 = DSINCAP)*DSINC INC) 

R23 = DCOS(AP)*DSINC INC) 

K33 = DCOSCINC) 

P ONDES GB ORIS Re 

9 = R21"I Т К22 1 БЕРЕ ЕК 

И ЕВЕ ВВ 

ДЕШ ЕМ 

END 


9222927227269225%921652259252 


DSINCLAN)*DCOS(ÀP) * DCOS(LAN)*DSIN(AP)*DCOS( INC) 
-DSINCLAN)*DSIN(AP) % DCOS(LAN)*DCOS(AP)*DCOS( INC) 
-DCOSCLAN)-DSTNQENS)) 
DSINCAPOSDSINOGI S) 

DCOS(AP) DSTA 
= DCOSC INC) 

R115 Р + К1270 + R13*W 

R2]*P + R220 ТЕО 

R3I*P т АЗОН 
URN 


uU n Hu n 


ROUTINE IJKPQWCLAN AP, ING ТЛ 


ге је У УУ У У) УУ У У У де У о и У те је у о У У је esee deese 


SUBROUTINE IJKRSW( LANAC по е ро ND 


THIS SUBROUTINE CHANGES FROM IJK COORDINATES TO RSW COORDINATES 


DOU 
ou 
R11 
R12 
I NS 
R21 
Кра 


КЕШ 
EAD 


BLE PRECISION 
R31,R32,R33,LAN,AL 
- DCOS(LAN)*DCOS(AL) - DSIN(LAN)*DCOS( INC)*DSIN(AL) 
DSIN(LAN)*DCOS(AL) * DSIN(AL)*DCOS(LAN)*DCOS( INC) 
DSIN(INC)*DSIN(AL) 

-DCOS( LAN)*DSIN( AL) -DSIN(LAN)*DCOS( INC)*DCOS( AL) 
-DSIN( LAN)**DSIN(AL) + DCOS( LAN)*DCOS( INC)*DCOS(AL) 
DSIN( INC)*DCOS( AL) 
DSIN( LAN) *DSIN( INC) 

-DCOS( LAN)*DSIN( INC) 

DCOS( INC) 

21191 + R12*J + R13*K 

R21*I + R22*J + R23*K 

R31*I 4 R32*J 4 R33*K 
URN 


Uu gg wu Ww H I 


ТАС, ви, ТВ БС MEOS 


0804440 
ОКВО44 50 
0КВ04460 
08804470 
ОКВО4480 
08804490 
08804500 
0КВ04510 
08804520 
08804530 
08804540 
08804550 
0804560 
ORB04570 
ORB04580 
0КВ04590 
08804600 
05804610 
0804620 
08804630 
ORB04640 
08804650 
0804660 
0КВ04670 
08804680 
0КВ04690 
ORB04700 
ОКВО4 710 
ОКВО4 720 
ORB04730 
ORB04740 
ORB04750 
ORB04760 
ОКВО47 70 
ОКВО4 780 
ОКВО4 790 
08804800 
ORB04810 
ORB04820 
ORB04830 
ОКВО4840 
0КВ04850 
ОКВ04860 
0КВ04870 
0КВ04880 
0КВ04890 
ORB04900 
0КВ04910 
ORB04920 
ORB04930 
ORB04940 
0КВ04950 
0КВ04960 
ОКВ04970 
ORB04980 


Veyezese УУ У И И И У У И У УУ УУ eee 


SUBROCPINE RSWIIKCLAN SAL, ING, R,S,W,1,J,K) 
d THIS SUBROUTINE CHANGES FROM RSW COORDINATES TO IJK COORDINATES 


Meera PRECISION MINC R S K I I KRIT RI2,R13 R21,R22,R23, 
+ R31,R32,R33,LAN,AL 


R11 - DCOSCLAN)*DCOS(AL) - DSINCLAN)*DCOSC INC )*DSIN(C AL) 
R21 = DSINCLAN)*DCOSCAL) + DSINCAL)*DCOS( LAN) *DCOSC INC ) 
R31 = DSINCINC)*DSINCAL) 

R12 = -DCOS(LAN)*DSIN(C AL) -DSIN( LAN) *DCOSC INC )*DCOS(C AL) 
R22 2 -DSINCLAN)*DSIN(AL) + DCOSCLAN)*DCOSC INC )*DCOS( AL) 
R32 = DSINCINC)*DCOS( AL) 

R13 = DSINCLAN)*DSINC INC) 

R23 = -DCOS(LAN)*DSIN(INC) 

R33 2 DCOS(INC) 


ПЕЕЛИ Уо ВІЗИ 
TEIR I RE К2255 + R23*W 
КИК Кк R32*5 + R33*W 
RETURN 

END 


Jesevesevesesesezesevevevedee deve vese eee eese de eee desee ye deve ve deze eyes eye ey eeveye deese ees eee ee eeesezeg e 


SUBROUTINE PQWRSW(TA,P,Q,W,R,S,WN) 
х THIS SUBROUTINE CHANGES FROM PQW COORDINATES TO RSW COORDINATES 


DOUBLE PRECISION P,Q,W,R,S,WN,R11,R12,R13,R21,R22,R23, 
+ — R31,R32,R33,TA 


Rll = DCOS( TA) 

[Ez DSIN(TA) 

КІЗ - 0.0 

R21 = -DSIN(TA) 

о = DCOS(TA) 

R23 = 0.0 

R31 = 0.0 

R32 = 0.0 

Roo = 1.0 

К = RII*P + R12*Q + R13*W 
бош  R21* РЕ. Во 
WN = R31*P + R32*Q +R33"W 
RETURN 

END 


УУУУ УСУ estes. se ә: - е: esee sese estos se 75: se ое УУУУ ЧОУ чоо очо а До Дора ада да а, сит ~“ 


SUBROUTINE RSWPQW(TA,R,S,W,P,Q,WN) 
* THIS SUBROUTINE CHANGES FROM RSW COORDINATES TO PQW COORDINATES 


О ЕВ ӘК SO Есиен КЛ ЕКІ RIS, R21,R22,R23, 
в вв. А 


R11 = DCOS(TA) 
R21 = DSIN(TA) 
R3i = 0.0 

R12 = -DSIN(TA) 


ORB04990 
ORBO5000 
ORBO5010 
ОКВ05020 
ORBO5030 
ORB05040 
ORB05050 
ORBO5060 
ОКВО5070 
ОКВ05080 
ОКВО5090 
ОКВО5100 
ОКВ05110 
ORBOS 120 
ОКВ05130 
ОКВО5140 
ORB05150 
ОКВ05160 
ОКВО5170 
ОКВО5180 
ORBO5190 
ORBO5200 
ОКВО5210 
ORB05220 
0КВ05230 
ОКВО5240 
ORB05250 
ORB05260 
ОКВО5270 
ORBO5280 
ORBO5290 
ORBO5300 
ОКВО5310 
ОКВО5320 
ОКВО5 330 
ORB05340 
ORBO5350 
ОКВ05360 
ОКВ05370 
ОКВ05380 
ОКВО5390 
ОКВ05400 
ОКВ05410 
ORB05420 
ORB05430 
ORB05440 
ОКВО5450 
ОКВО5 460 
ОКВ05470 
ОКВ05480 
0КВО5490 
ОКВО5500 
ОКВ05510 
ORBO5520 
ORBO5530 
ОКВО5540 


R22 - DCOS(TA) 

R32 = 0.0 

R13 = 0.0 

Rese = 0 9 

Ree = lee 

ра КА К ERI 
Q = R21*R 4852255 4752354 
WN = R31*R + R32*S +R33*W 
RETURN 

END 
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TORE ELEMENTS IN ARRAYS 
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SUBROUTINE STORE(RI,RJ,RK,R,TA,RIRAY ,RJRAY , RKRAY , RARAY , TARAY ,„МИМ, 
1.,AP,AINRAY , ARRAY, TE ПЕКАТ! 

THIS SUBROUTINE STORES THE POSITION AND ELEMENTS SIN AKRAYS sit 

SINGLE PRECISION СӘКІ БОВ ВОН ING 


DOUBLE PRECISION Rigko, RRVRG TA СТАРТ 


DIMENSION RIRAY( 500) ,RJRAY(500) ,RKRAY( 500) , RARAY( 500) , TARAY( 500), 


AINRAY( 500) ,APRAY( 500) , TIMRAY( 500) 
RIRAYC NUM) = SNGLCRI) 
RJRAY(NUM) = SNGL(RJ) 
RKRAYCNUM) = SNGLCRK) 
RARAY(NUM) = SNGLC(R) 
TARAY(NUM) = SNGL(TÀ) 


AINRAY(NUM) 9» SNGL(I) 

APRAY( NUM) = SNGLCAP) 

TIMRAY(NUM) = SNGLCTT) 
RETURN 

END 


ye vevevezesevezevedeves ese ese esee deve ye Je vesc Jesedes eye ede de esee deve de deve ye des eye deese ye eese esee 


INITIAL POSITION, VELOCITY 


yeyesese ev eese esee desee evedese dee dee eeevyedes e yevey ez eese ev esee se cesesedevcsesesicseses Хе о О ЛГУ: КЛ УЛ 


se 


SUBROUTINE INPUTSCRI,.RIJISRKOR VIVUS КЕЕ, 

THIS SUBROUTINE GIVES IHE USER ATC FOT ВЕСЕ ТЕК THe 
INITIAL POSITION AND VELOCITY VECTOR OR TO LET THE PROGRAM 
CALCULATE THE INITIAL POSITION AND VELOCITY FROM USER PROMPTED 


INPUTS 

SUBROUTINES CALLED FROM THIS SUBROUTINE: 

INELTS = Prompts USER FOR ORBITAL ELEMENTS 

IPOS = PROMPTS USER FOR INITIAL POSITION (IJK) 
IVEL = PROMPTS USER FOR INITIAL Velocity (IJK) 


DOUBLE PRECISION КК, ЕКТ 
CHARACTER* 1, QUIT 


PROJET USER FOR METHOD: TO ЕМЕ ИРЕ 
PRINT, IN WHICH MANNER WOULD YOU LIKE TO INPUT THE INITIAL’ 


(„2 
СМ 


ОКВО5550 
ОКВО5560 
0805570 
ORB05580 
0КВ05590 
ОКВО5600 
ОКВО5610 
ORB05620 
ОКВО5630 
ОКВО5640 
ОКВО5650 
ОКВО5660 
ORB05670 
ORB05680 
0КВ05690 
ОКВО5 700 
ОКВО5 710 
ОКВО5 720 
ОКВО5 730 
ORB05740 
ORB05750 
ORB05760 
ORB05770 
ОКВО5 780 
ОКВО5 790 
ОКВО5800 
ORBO05810 
ОКВО5820 
ОКВ05830 
ОКВО5840 
ORBO5850 
ORBO5860 
ОКВО5870 
0КВ05880 
ОКВО5890 
ОКВО5900 
ОКВО5910 
ORBO5920 
ОКВО5930 
ОКВО5940 
0КВ05950 
ОКВО5960 
ORB05970 
ОКВО5980 
ОКВО5990 
ОКВ06000 
ORB06010 
ORB06020 
ОКВО6ОЗО 
ОКВО6О040 
ОКВ06050 
ОКВ06060 
ОКВ06070 
ОКВ06080 
ORB06090 
ORBO6100 


4... 


180 


ыа ае а eg 
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PRINT*,'POSITION AND VELOCiTY OF THE SATELLITE?' 


PRINT," 1: BY Inputting THE INITIAL POSITION AND VELOCITY' 
PRIST, VECTORS IN THE PERIFOCAL COORDINATE SYSTEM (IJK)' 
БЕР 2. BY LETTING THE SATELLITE BE PLACED ON THE "I'"' 
PRINT*,' ANTSOOENENE CIOJERSISYSTENMCAT A DESIRED RADIUS OF' 
PRINT*,' PERIGEE(RP) AND INPUTTING EITHER A DESIRED RADIUS' 
PRINTS OF APOGEE(RA), A DESIRED ECCENTRICITY(E), OR THE' 
PRINT*, ' DESIREU) VELOCI UY AW THAT RADIUS. AND A DESIRED 
PRINT*, | PACLINATICNG 

PRINT*,' SEDI TT: 

PRINT РЕМ ОЛЛО eo ' 

READ* , ICHC 


PRINT*, ICHC 
CALL EXCMS( CLRSCRN') 


USER INPUTS POSITION AND VELOCITY VECTORS 
DE TCRO TEOT) THEN 

ФАП POstni RI] RKR) 

АБУ УТ, УЛ, УК. У. в, М0) 


АВВ ТЕСТО ORBITAL ELEMENTS ТО GET POSITION AND VELOCITY 
ВЕ ВОН ВО 2) THEN 
CrbleiNELISCRI,RIZRK,R,VI,VJ,VK,V,MU;P1) 


STOP PROGRAM 

Meee CTCHO .EQ. 3) THEN 
ОТ = № 

ENSE 
PRINT*,'INVALID ENTRY! 
COTO 195 

ENDIF 

RETURN 

END 


TRY AGAIN! ' 


Sees were wes es eee eee eee eve eee e eee eve eve ee eve eee 
SUBROUTINE IPOSCRI,RJ.RK,R) 

THIS SUBROUTINE ASKS THE USER FOR THE INITIAL POSITION OF THE 
Satellite IN GEOCENTRIC-EQUATORIAL COORDINATE SYSTEM 

IunsHsE PRECISION RI,RJ,RK,R 
CHARACTER*1, CHOICE 

LOGICAL CORREC 

CORREC = .FALSE. 


PROMPT USER FOR VELOCITY VECTOR 
IF(. NOT. CORREC) THEN 
CAME ence CELRSCRY > 
PRINT*,'ENTER RADIUS VECTOR VALUES IN КМ 
PRINT*,'RADIUS OF THE EARTH = 6400 ҚКМ! 
CORREC = TRUE. 
PRINT*,'ENTER RI :' 
READ* ,RI 
PRINT“, RI = ОРТ, КМ. 
PRINT*,' ENTER RJ :' 


ОКВ06110 
ОКВ06120 
08806130 
ORB06140 
ОКВ06150 
ОКВ06160 
ОКВ06170 
ОКВ06180 
ОКВ06190 
ОКВ06200 
ORB06210 
ORB06220 
ОКВО6230 
ORB06240 
ORB06250 
ОКВО6260 
ОКВО6270 
ОКВО6280 
ОКВО6290 
ОКВ06300 
ОКВ06310 
ОКВ06320 
ОКВ06330 
ОКВ06320 
ORB06350 
ОКВ06360 
ОКВ06370 
ОКВО6ЗВО 
ОКВ06390 
ORBO6400 
ORBO06410 
ORB06420 
ORB06430 
ОКВО6440 
ОКВО64550 
ОКВО6460 
ОКВО64 70 
ОКВ06480 
ОКВО6490 
ОКВ06500 
ОКВ06510 
ОКВ06520 
ОКВ06530 
ОКВ06540 
ОКВ06550 
ОКВ06560 
ОКВО6570 
0КВ06580 
ОКВ06590 
ОКВ06600 
ОКВ06610 
ОКВ06620 
ОКВ06630 
ORBO6640 
ORB06650 
ОКВ06660 


вр >ò 9 e 
m oo 


IU 


READ* ,RJ 
PRINTS C RI SARI рай 
РЕШЕ БУРЕК Бенин 
КЕАО „ЕК 

PRINT, ЕК = ЕКОШ 


CALCULATE TOTAL R 

в = DSQRT( (RI***2 ) + ( RJ**2) + CRK**2) ) 

PRINT*,'R = °,R,' KM 

ЈЕ (Е БЕ. 640020 ОРЫ 
PRINT*, RADIUS TO SMALL! ! 
GOTO 180 

ENDIF 


ENTER NEW VALUES!! ' 


CHECK WITH USER THAT Values ARE CORRECT 
PRINT*,'ARE THESE VALUES CORRECT?' 


PRINT ENTER "Y OREN: 
READ* , CHOICE 
CHONG hE =, 


PRINT* , CHOICE 

LF (CHOICES EO ms) sre: 
CORREC = . TRUE. 

ENDIF 

GOTO 180 

EADIE 
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SUBROUTINE ПЕ Ув.) 
THIS SUBROUTINE ASKS THE USER FOR THE INITIAL VELOCIT: Obamas 
Satellite 


DOUBLE PRECISION VI,VJ,VK,V,R,VCIR, VMAX MU 


CHARACTER*1, CHOICE 
LOGICAL CORREC 
СОККЕС = .РАШЗЕ. 


CALCULATE ESCAPE VELOCITY AND CIRCULAR VELOCITY AND PROMPT USER 
FOR VELOCITY VECTOR 
IF(. NOT. CORREC) THEN 

CALE EXENSC CLRSCRN 

VCIR = DSQRT(MU/R) 

VMAX = DSQRT( (2. 0**MU)/R) 

PRINT*, CIRCULAR VELOCITY = OIR KM/SEC 

PRINT*, ‘MAXIMUM VELOCITY = ', VMAX, 'KM/SEC' 

CORREC = .TRUE. 

PRINT*,' ENTER VELOCITY VECTOR IN (KM/SEC)' 


PRINT“, 'ENTER VI : ' 


READ* , VI 

PRINT Ul = ДЕНЕ 
PRING ТЕКТЕН О 
READ , VJ 


2 
ur 


ОКВО66 70 
ОКВО66ВО 
ОКВО6690 
ОКВО6 700 
ОКВО6 710 
ОКВО6 720 
ОКВО6 730 
ОКВО6740 
ОКВО6 750 
СКВ06760 
ОКВО67 70 
ORB06780 
ORB06790 
ORBO6800 
0КВ06810 
ORB06820 
ORB06830 
ORB06840 
0КВ06850 
ORBO06860 
0КВ06870 
0КВ06880 
0КВ06890 
0КВ06900 
ОКВ06910 
0КВ06920 
ORBO06930 
CRB06940 
0КВ06950 
0ОКВ06960 
0КВ06970 
ORBO06980 
0КВ06990 
ОКВ07000 
ОКВ07010 
ORBO7020 
ОКВ07030 
ОКВ07040 
ORBO7050 
ОКВО 7060 
ORBO7070 
ОКВО 7080 
ORBO7090 
ORBO7100 
ORBO7110 
ORBO7120 
ORBO7130 
ORBO7140 
ORBO7150 
ORBO7160 
ОКВО7170 
ОКВО7180 
ORBO7190 
ORBO7200 
ORBO7210 
ORBO7220 


PRINT*,'VJ - ',VJ,'KM/SEC' ОКВО 7230 


PRINT, ‘ENTER VK : | ORBO7240 

READ* , VK ORB07250 

Hune у = VK. NM/SEC- ORBO7260 

ORBO7270 

* CABCULATE TOTAL VELOCITY (v) ORB0O7280 
V - DSQRT((VI|*2) 9 (VJ9*2) - (VK**2)) ORB07290 
PRINT*,'V - ',V,'KM/SEC' ORBO7300 

ORBO7310 

* CHECK WITH USER THAT VALUES ARE CORRECTS ORBO7320 
PRINT*,'ARE THESE VALUES CORRECT? ' ORB07330 

PRINT*,' ENTER "Y" oR "Ww" ;' ORB07340 

READ* , CHOICE ORB07350 

CRO CNS. У. ORB07360 
PRINT*,CHOICE ORB07370 

IF (CHOICE. EQ. ' Y' ) THEN ORBO7380 

СОККЕС - .ТЕСЕ. ORB07390 

ENDIF ORBO7400 

IF (V .GE. VMAX) THEN 0КВО7410 
PRINT*, ' VELOCITY IS GREATER THAN THE ESCAPE VELOCITY!!' ORB07420 

PRINT*,' RE-ENTER VELOCITY!!! ' ORB07430 

CORREC = . FALSE. ORBO7440 

ENDIF ORB07450 

GOTO 190 ORBO7460 

ENDIF ORBO7470 
RETURN ORBO7480 

END 0ВВ07490 
0ЕВО7500 

Уа УУ У ТУ У УЖИ ПУК и хе пху ORBO7510 
ORBO7520 

ENGROUTINE INELBTSCRI, RJ, RK, R, VI, VJ, VK, V,MU,;PI) ORBO7530 

* SATELLITE PLACED ON 'I' AXIS AND USER SUPPLY ORBITAL ELEMENTS TO ORBO7540 
* Gael. 1TIAL POSITION AND VELOCITY ORB07550 
ORBO7560 

DOUBLE PRECISION RI,RJ,RK,R,VI,VJ,VK,V, MU, I,ENR,A,E,RP,RA,PI,VMAX ORBO7570 
CHARACTER*1,CHOICE ORB07580 
ОЕВО7590 

* PROMPT USER FOR PERIGEE RADIUS ORBO7600 
198 | PRINT*,'ENTER RADIUS OF PERIGEE(RP) IN (KM), FOR EXAMPLE: ' ORBO7610 
PRINT*,'LOW EARTH ORBIT (LEO), RP = 6600.0 KM' ORBO7620 
PRINT*,'GEOSYNCROUNOUS ORBIT, RP - 42241. 1 КМ! ОЕВО7630 
PRINT*,'ENTER RP: ' ORB07640 
PRINT*,'"RP" MUST BE > 6400KM' ORB07650 

READ* , RP ORBO7660 
PRINT*,RP ORB07670 
ORB0O7680 

* CHECK FOR VALID RADIUS 0ВВ07690 
IF (RP .LT. 6400.0) THEN ОВВ07700 
PRINT, 'YOUR "ЕР" IS TO SMALL!! ' ORBO7710 

GOTO 198 ORBO7720 

ENDIF 0ВВ07730 
ORBO7740 

х PROMPT USER FOR TYPE OF INPUT ORB07750 
PRINT, DO YOU WANT TO ENTER THE ECCENTRICITY (E), ' ОВВ07760 
PRINT*, RADIUS OF APOGEE (RA), OR VELOCITY (V)?' ORBO7770 


Е EMER ER’. ORV: ORBO7780 


(,2 
MO 


READ* ,CHOIC 
PRINT* “CHOICE 
CALL EXC {S('CLRSCRN' ) 
USER ENTERS Eccentricity AND SEMI-MAJOR AXIS, ENERGY AND VELOCITY 
IS CALCULATED IN THAT ORDER 
IF (CHOICE .EQ. 'E') THEN 

PRINT*,'ENTER ECCENTRICITY (E): ' 

РЕТМТ“, 0.0 <= E < 1.0' 

READ* ,E 

PRINT*,E 


* CHECK FOR VALID ECCENTRICITY 
IF CCE .LT. 0.0) ОЕ (Е о оловни 
PRINT*,'INVALID "E'' 
GOTO 198 


в 
А = ДИ E) 

ENR = М0/(2. О АД) 

ү = DSORT( 2 “CENR+(MU/RP) )) 


5 USER INPUTS RADIUS OF APOGEE AND ECCENTRICITY IS CALCULATED 
хе THEN SEMI-MAJOK AXIS, ENERGY AND THEN) VELOCITY: 
БӘЗЛЕР (СНОЛОНВ ЕП БЕТІНІН 
PRINT*, ‘ENTER RADIUS OF APOGEE (RA) IN K™: | 


PRINT*,'"RA" MUST ВЕ > БР, "ЕР < " „БР 
READ* ,RA 
PRINT*,RA 


CHECK FOR VALID RADIUS OF APOGEE 
TEC Avec TUE RP) THEN 
PRINT, YOUR” RA’ IS TOSI LINE 
GOTO 198 
ENDIF 
(RA-RP)/(RA+RP) 
ОС) 
А) 
DSQRT(2*(ENR*(MU/RP))) 


<-> uit 
А 
| и | 


> USER INPUTS MAGNITUDE OF VELOCITY, PROGRAM PROVIDES CIRCULAR 
* AND ESCAPE VELOCITY FOR COMPARISON AND TO CHECK FOR VALID 
vs INPUTS 
ELSEIF (CHOICE .EQ. 'V') THEN 
PRINT**, 'ENTER VELOCITY IN KM/SEC: ' 
PRINT*,'THE MINIMUM VELOCITY ALLOWED IS FOR A CIRCULAR ORBIT’ 
VCIRC = SQRT(SNGL(MU/RP)) 
PRINT*, ORBIT. V€Cireular) = eel RG aro ou 
VMAX = DSQRT(2**(MU/RP)) 
PRINT*,' THE MAXIMUM VELOCITY < ',VMAX,' KM/S' 
READ* ,V 
PRINT*,V 
IF (V . LT. VCIRC) THEN 
PRINT*,'VELOCITY TO SMALL!' 
GOTO 198 
ENDIF 
F (V .GE. VMAX) THEN 


ОКВО7 790 
ORBO7800 
ORB07810 
ОКВО 7820 
ОКВ07830 
ORB07840 
ОКВО7850 
ОКВО 7860 
ORBO7870 
ORBO7 880 
0КВ07890 
ОКВО 7900 
ORBO7910 
ORB07920 
ORBO7930 
ORB07940 
0КВ07950 
0КВ07960 
0КВ07970 
ORBO7980 
ОКВО 7990 
0КВ08000 
ОКВ08010 
ORB08020 
ORB08030 
ORB08040 
ORB08050 
ОКВ08060 
ORB08070 
ОКВОВОВО 
ORBO8090 
0КВ08100 
ОКВ08110 
0КВ08120 
ОКВ08130 
ORB08140 
ORB08150 
ОКВ08160 
ORB08170 
ORB08180 
ORB08190 
ORB08200 
ORB08210 
ORB08220 
ORB08230 
ORB08240 
0КВ08250 
0КВ08260 
ОКВОВ2 70 
ORB08280 
ORB08290 
ORB08300 
ORB08310 
ORB08320 
0КВ08330 
ORB08340 


+ 


РЕТХТ“,' УЕГОСТТУ ТО СКЕАТ!!' 
GOTO 198 
ENDIF 
ELSE 
PRINT^,'INVALID ENTRY! TRY AGAIN' 
GOTO 198 
ENDIF 


INCLINATION NEEDED TO GIVE Velocity A Direction 
PRINT*, ‘ENTER INCLINATION (I) IN DEGREES: ' 
READ* , I 

PRINT*,I 

ШЕСГЕТ/1280. 071 

VK = V*DSINCI) 

VJ 2 V*DCOS(I) 

VI 0. 0 


ВЕВТОЬ VECTOR SET 


д) 
Y C 
ШІ 
© 
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ee OUTINE CALCEL(CRI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E,A,1I,LAN, 


THIS SUBROUTINE CALLS THE INDIVIDUAL SUBROUTINES TO CALCULATE THE 


ВЕ Е РД, НН, НМ 


ORBITAL ELEMENTS 


THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES(RETURNED VALUES) 


ENERGY = ENERGY PER MASS (ENR) 

ANGMOM = ANGULAR MOMENTUM (H,HI,HJ,Hk) 
NODE = NODE VECTOR (N,NI,NJ,NK) 

LATREC = SEMI-LATUS RECTUS (P) 

ECC = ECCENTRICITY (E,EI,EJ,EK) 

SMAXIS = SEMI-MAJOR AXIS (A) 

ЕСІ = INCLINATION (I) 

ASNODE = LONGITUDE OF ASCENDING NODE (LAN) 
ARP = ARGUMENT OF PERIGEE (AP) 

IJKPQW = 'IJK' SYSTEM TO 'PQW' SYSTEM 
TANOM = TRUE ANOMALY (TA) 

ARLAT = ARGUMENT OF LATITUDE (AL) 

LONPER = LONGITUDE OF Perigee (LP) 

TLON = TRUE LONGITUDE (TL) 

PERIOD = PERIOD (PER) 

bees. = ECCENTRIC ANOMALY (EA) 

MEANMO = MEAN MOTION (MM) 

MEANAN = MEAN ANOMALY (MA) 

TFLGHT = TIME OF FLIGHT (TF) 

В ВЕБ ТО В. БЕ РТ В Е АТ ПАК, АТ, 


4- 


ORB08350 
ORB08360 
ORB08370 
ORB08380 
ORB08390 
ORB08400 
ORB08410 
ORB08420 
ORB08430 
ORB0O8440 
ORB08450 
ORB08460 
ORB08470 
ОКВОВАВО 
05808490 
ORBO8500 
ОКВО85 10 
ОКВОВ5 20 
0КВ08530 
ORB08540 
ORBO08550 
ОКВО8560 
ОКВОВ5 70 
0КВ08580 
0КВ08590 
ORB08600 
ОКВ08610 
ORB08620 
ОКВОВ6ЗО 
OR308640 
ORB08650 
ORBO8660 
ORB08670 
ОКВ08680 
0КВ08690 
ORB08700 
ОКВОВ7 10 
ОКВОВ 720 
ОКВ08730 
ОКВ08740 
ОКВОВ 750 
ОКВ08760 
ОКВ08770 
ОКВ08780 
0КВ08790 
0КВ08800 
0КВ08810 
0КВ08820 
ORB08830 
ORB08840 
ORB08850 
ORB08860 
ORB08870 
ОКВОВВВО 
0КВ08890 
ORBOS900 


' LU 
-- 5-759 


* % . . т % * = , . е As [J * % % *. + * . ? 
еее сеет 


Р,ТА,РЕК,ЕА, МА. АР, ТР. АТН, ВЕ, Но пика, 
Т.В. ВО, КУ, ор. О 


САБ ЕХЕБОТ(( А К ЧО ENRI 
САБО О Ve oe 
CALL NODECHI,HJ, MINAS SESS 

CADESEATREGCHAP SMS 

CALL ECCCRI ЕТ ККр МП УТО МЕК V.EDSEJSOERS SEDED 
CALL SMAXISCMU,ENR, A) 

САТЕ ОИСЕ Н SN 


SPECIAL CASE IF INCLINATION = 0.0 
IF (ТО ОО ОЕ 

CALL ASNODE(NI,N,LAN,NJ,PI) 

CALL ARP(NI,NJ,N,EI,EJ,EK,E,AP,PI,NP,NQ, LAN) 
ELSE 

LAN = 0.0 

AP = 0.0 
ENDIF 


COORDINATE TRANSFORMATION OF 'R' AND 'V' VECTORS 

CALL IJKPQW(LAN,AP,I,RI,RJ,RK,RP,RQ,RW) 

CALL IJKPQW(LAN,AP,I,NI,NJ,NK,NP,NQ,NW) 

CALL TANOMCEILEJ EK E RD ORI RR REA ROW R QUE UK ASI 


SPECIAL CASE FOR Inclination < 0.0 
LR QUE NE. ЖО ЈАНЕ А 
CALESERDATGNINJONSENSSRITRJORA RADI AA 
ВЕ 
AL = TA 
EN Uae 
CALL LONPERC LAN AP VLE) 
CALL TLONCLAN, AP ,TA,TL) 
CABDSPERIODSASBEDMPDI МО) 
CALL ECCAN(E, TA EA PI) 
Сабо Еа DTE 
CALL MEANAN(EA,E,MA) 
CALL TEBGHTCOHHISSASTES 
RETURN 
END 


o a an n n ona nl a n'a nl n an a alata йо 5 ма та Ја Ја а о зъба ааа? 
AN 6*8 48 Фу зъ Фу Фу съ 272: PN FDEN АА КОТА И ХУ ААА ИЯ 


PPP PT Pm 
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SUBROUTINE ENERGY(V,R,MU,ENR) 
THIS SUBROUTINE CALCULATES THE ENERGY OF THE ORBIT 


DOUBLE PRECISION V,R,MU,ENR 
ЕКСЕ ЕН 
RETORN 

EAD 


се узео у О У ет  е 


SUBROUTINE АХСЧОМСЕТ. КЛ КА МП РР. ДВС Туле ОНА И) 
THIS SUsROUTINE CALCULATES THE ANGULAR MOMENTUM 


ЕЎ БЯ 


08808910 
ORB08920 
ORB08930 
ORB08940 
0КВ08950 
0КВ08960 
ОКВО8970 
ОКВО8980 
ОКВО8990 
ОКВО9000 
ОКВ09010 
ОКВ09020 
ОКВ09030 
ОКВО9040 
ORB09050 
ОКВО9060 
ОКВО9070 
ОКВ09080 
ОКВО9090 
ОКВО9100 
ОКВО9110 
ORB09120 
ORB09 130 
ОКВО9140 
ОКВО9150 
ORBO9160 
ОКВ09170 
ОКВ09180 
ОКВО9190 
ОКВО9200 
ORB09210 
ORBO9220 
ORB09230 
ORB09240 
ОКВО9250 
ОКВО9 260 
ORB09270 
ОКВО9280 
ORB09290 
ОКВ09300 
ОКВ09310 
ORB09320 
ORBO09330 
ORB09340 
ОКВО9350 
ORB09360 
ORB09370 
ОКВО9380 
ORB09390 
ORB09400 
ORB09410 
ORB09420 
ORB09430 
ORB09440 
ORB09450 
ОКВО9460 


Jegede yes 
СУТЕК 


Dever PRECISION KI,RJ;RK,VI,VJ,VK,HI,HJ.HK,H 


ОЕК) (ВК * VJ) 
о) + (ВР УМЕО 


HK О ТГ) 
EEREODSORTCCHIS*2) + (HI**2) + (HR **2)) 
RETURN 
END 
Joevedesesevez eyes edes esevyesevesevevezes eve eee desee eye ee eyezevevevyedese eye eveve see veyevfevevese sev eses eve vez Ж 


SUBROUTINE NODE(HI,HJ,NI,NJ,NK,N) 
THIS SUBROUTINE CALCULATES THE NODE VECTOR 


BOURLE PRECISION HI,HJ,NI,NJ,NK,N 


222 
e 
нии 
күз 
+4 


КК 0.0 

DE-SDSORTCCNIS*2) - (NJS*2)) 
RETURN 

END 


ИУ УУ И ИУ УУ У У У У У У У У У У У У УУ УУ: 


SUBROUTINE LATREC(H,P,MU) 
Mais SUBROUTINE CALCULATES THE SEMI-LATUS RECTUM 


DOUBLE PRECISION H,P,NU 
В НО 


RETURN 
END 


ene ae lon! m» mm “--. а о = a's a's “.-....... А ъв а е ае во а ъв във в ъв ъбо аа е во ааъ о ъв ъйо ь! ж ъ' = йо ъв ъ е аа ъйо ъс ъй в ъс >: “..... we m ~“. 
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ДИ ОППВ ВОС КТ КЛ ККК. МЕ МУ, МК М БТ ЕЛ ЕК Е МО) 
THIS SUBROUTINE CALCULATES THE ECCENTRICITY 


Ша ше PRECISION RI, RI, RRR, Vigvs, VK,V,B1,EJ,Ek,£,MU,DOT 


CALCULATE DOT PRODUCT OF 'R' AND 'V' VECTORS 
DOT = (RI*VI) + (RJ*VJ) + (RK*VK) 
= Or OO Uae (CCV т (МШК) ЕТ = (рот) 
ЕКОО ОО НО) (ОСУ На) = CHU/R)) ВЛ - (ОГУЗ 
= (1.0D-00/MU) * (((V**2) - (MU/R)) * RK - (DOT)*VK) 


RETURN 
END 


ee E E АЕ а А Рана EEA TI AS TUA T T A O A TA а 


ОДРА ОВ ТОМЕ 5.18 ТОМО, ЕЗЕЈА) 
THIS SUBROUTINE Calculates THE SEMI-MAJOR AXIS 


ORBO9470 
ORB09480 
ORB09490 
ORB09500 
ORBOOS 10 
ORB09520 
ORBO2530 
ORB09540 
ORB09550 
ОКВО9560 
ОКВО9570 
ORB09580 
ОКВО9590 
ОКВО9600 
ORB09610 
ORB09620 
ORB09630 
ОКВО9640 
ORB09650 
ОКВ09660 
ОКВО9670 
ORB09680 
ORB09690 
ОКВО9 700 
ORBOS TIO 
974521012, 72720) 
066099790 
ОКВО9740 
ОКВО9 750 
ОКВО9760 
ОКВО9 7 70 
ОКВО9 780 
ОКВО9 790 
ОКВ09800 
ORB09810 
ORB09820 
ORB09830 
ORB09840 
ORB09850 
ОКВО9860 
ОКВО9870 
0КВ09880 
ORB09890 
ORB09900 
0КВ09910 
0КВ09920 
0КВ09930 
ОКВО9940 
ОКВО9950 
ORB09960 
ORB09970 
ORB09980 
00509990 
ОКВ10000 
ORB10010 
ORB10020 


DOUBLE PRECISION MU,ENR,A 


-— VETE атқ” 
А = = О а ЕН) 
уз Рут + ы 
RETURN 
END 
УУ те Akr k У У СК тес ее еее еее еее ет еее ЗҮ ТСС 


977929622229: 


SUBROUTINE INCLCHK,H,I,PI) 
THIS SUBROUTINE CALCULATES THE INCLINATION 
'I' ALWAYS LESS THAN 180 DEGREES 


DOUBLE PRECISION ЕТ 


I = DACOS( HK/H) 
RETURN 
END 


І29:2:9:92222%2692992965922:92929396769:596269:32569с6969:369с9769925269252925225797 


SUBROUTINE ASNODE(NI ,N,LAN,NJ,PI) 
THIS SUBROUTINE CALCULATES THE LONGITUDE OF THE ASCENDING NODE 
[END £26 THEN 'LAN' < 180 DEGREES 


DOUBLE PRECISION NI,N,LAN,NJ, PI 


LAN = DATAN2(NJ,NI) 

ТЕ CLAN LI. 90.0) ЕК 
LAN ="(2P1)- + GAN 

EXDIC 

RETURN 

EU 


“ааа! а ва ав ъв ива ве! а аа а ааа ба би? в ъв а в ъв alan aa as aula uan aa ue ae о агага ъа ога ъа њо њо иаа оо аба ъбл а ва ъв чь! @ ый ай ван tanta ela 272222... ae 
и пити 


SUEROUTINE АКРОМЕМЈ МЕТ ЕЈ,ЕКЈЕ APTER LAN) 
THIS SUBROUTINE CALCULATES THE ARGUMENT OF Perigee 
IF 'EX' GREATER THAN O THEN 'AP' « 180 

VARIABLE TEMP USED AS A Temporary VALUE FOR ARCTAN 


DOUBLE PRECISION NI,NJ,N,EI,EJ,EK,E,AP,PI,NQ,NP, TEMP, LAN 


IF (CED ВО О И Ея 
AP = 0.0 
ELSE 
TEMP = DATAN2(EJ,ET) 
ТЕ (TEMESI TELAN) THEN 
AP = TEMP - LAN 


(EJ EQ: TOT ODDA THEN 


> 
25) 
|| 


РЕ) они А нем БОЦЕ) 
ENDIF 
Е ар. 0 ODE 
AP = (ӘЖЕМЕ р 
Ваше 
TF (АР СТ. Сара НАЈ 
Ар т) 


ОКВ10030 
ORB10040 
ORB10050 
ORB10060 
ОКВ10070 
ОКВ10080 
ОКВ10090 
ОКВ10100 
ORB10110 
ORB10120 
ORB10130 
ORB10140 
ORB10150 
ORB10160 
ORB10170 
ORB10180 
ORB10190 
ORB10200 
ORB10210 
ORB10220 
ORB10230 
ORB10240 
ORB10250 
ORB10260 
ORB10270 
ORB10280 
ORB10290 
ORB10300 
ORB10310 
ORB10320 
ORB10330 
ORB10340 
ORB10350 
ORB10360 
ORB10370 
ORB10380 
ORB10390 
ORB10400 
ORB10410 
ORB10420 
ORB10430 
ORB10440 
ORB10450 
ORB10460 
ORB10470 
ORB10480 
ORB10490 
ORB10500 
ORB10510 
ORB10520 
ORB10530 
ORB10540 
ORB10550 
ORB10560 
ORB10570 
ORB10580 


BANDIT RE i 

БЕШЕ =e i 

ОКА 351 

END nel 

ШЕ 22205792225522227:5:72227257%22 TUI стр vero rdiet E IS аА = 231 

КВ] 

SUBROUTINE TANOM(EIL,EJ,EX,E,RI,RJ,RE,RP,RQ,REW,R, VI,VJ,VK, 281 

+ ТЕ,РІ RET 

F THIS SUBROUTINE CALCULATE TRE “PAGE” enoma ly RS 

E ~ ОО М) > 0 Тлах Је < 160 DEGREE 581 

БВ 1 

Ши ВЕ РЕБОТУТОХ“ DOT, EI, EJ, EK, E RI RI RER: VI У ЈЕ ДА PT, R81 

+ — RP,RQ,RW 231 

R31 

TA = DATAN2(RQ,RP) КВт 

ШЕ (ТА .17.0.0 ) ТИЕК РБ! 

1 Cee је RET 

ЕЕ RE! 

RETURN RB! 

END КЕ] 

=D i 

о о Не ета а асса аа 281 
В ОВИМЕ ARLATI(MI.NJ,NA,N,RI,RJ,RASS,AL,PI,TÀ,A2) z 
С ВЫ SUBROUTINE B. G THE ARGUMENT CF LATETUDE x 

P ШО > 0) TRENPAE < 150 DEGREES 
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SEBROUTINE ТОМРЕК( ТАМ АР ,„ТР) 
E BUS SUBROUTPNE CALCULATES THE LONCITUSE CF PERIGEE 
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DOUBLE PRECISION LAN ŠP, LP 
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SUBROUTINE TLONCLAN,AP,TA,TS 
= Mes SUBROUTINE CALCULSIES GRE IxUE LONGLIVDE ай EPOCI 
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DOUBLE PRECISION пак ще, РА, тт, 
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TL = AP + LAN : 
ТЕК О5Е11920 
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SESEOUTINE РЕРТОО А РЕК УР О) 
x THIS SUERGUTING CALCU ШЕСТЕ БЕРИ 


DOUBLE PRECISION ASBESSEISTIU 

PER = 2. 0D+00%*( PI)*DSQRT(( A***3) /MU) 
RETURN 

END 
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SUBROUTINE ЕССАМСЕ,ТА,ЕА РГ) 
* THIS SUBROUTINE CALCULATES THE ECCENTRIC Anomaly 


DOUBLE PRECISION E,TA,EA,PI 
EA = DACOS((E + DCOS(TA))/(1. OD+00 + E*DCOS(TA))) 


ТЕЧ ОТА СТ. РВ 
БА = ГГ 


у ит у Ло од Хо У 


SUBROUTINE МЕАММО(А,ММ,М0) 
THIS SUBROUTINE CALCULATES TRE. MEAN OT TON 


DOÜDSLE PRECISION Agta 
= DSORT MU Ш) 


RETURN 
END 
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SUBROUTINE MEANANCEA,E,MA) 
THIS SUBROUTINE CALCULATES THE MEAN Anomaly 


DOUBLE PRECISION BALE MA 
MA = EA - E*DSIN(EA) 
RETURN 

END 


Зее еее еее еее ее tededededevesicsc reve ve Fock ses ои dede reve se seve te sede reve 


SUBROUTINE TFLGHT(MM,MA,TF) 
т: THIS SUÜBROUTINE CABGBEATES THE ГІМЕ СЕЗРЕТОЛІ 


DOUBLE PRECISION MM,MA,TF 


ТЕ А 


4^ 


ORB11140 
ORB11150 
ORB11160 
ORB11170 
ORB11180 
ORB11190 
ORB11200 
ORB11210 
ORB11220 
ОКВ11230 
ОКВ11240 
ОКВ11250 
ORB11260 
ORB11270 
ORB11280 
ORB11290 
ORB11300 
ORB11310 
ORB T1320 
ORB11330 
ORB11340 
ORB11350 
ORB11360 
ORB11370 
ORB11380 
ORB11390 
ORB11400 
ORB11410 
ORB11420 
ORB11430 
ORB11440 
ORB11450 
ORB11460 
ORB11470 
ORB11480 
ORB11490 
ORB11500 
ORB11510 
ORB11520 
ORB11530 
ORB11540 
ORB11550 
ORB11560 
ORB11570 
ORB11580 
ORB11590 
ORB11600 
ORB11610 
ORB11620 
ORB11630 
ORB11640 
ORB11650 
ORB11660 
ORB11670 
ORB11680 
ORB11690 


RETORN ORBII700 

EXD ORBA #719 

ORE 39720 

ППИ ит РТ ТУР ТЕТ У И РУТЕ У ТЕ УТО ЕУ ТЕТ ТЕ И ТЕ ЕТ ТТ УУ ИЕ ОУ И ИЕ ти ЈДЕ21 1730 
4 BALCULATE UNPERIURBED СВЕТЕ ORBA 1740 
2595975555579: суви титули пп пепипи 05811750 
05811760 

SUBROUTINE UNPRET(DT,PER,AL, LAN, AP,1,R1,2J,2X,R, 25311770 
БЕ vr We V, MU. bL E A E N TA P MINE EN, 08811780 

E i Ts NUM,RIRAY,RJRAY,RKRAY,RARAY,TARAY,AINRAY,APRAY,TIMRAY, 05812790 

T TT) 05811800 

d THiS SUEROUTINE CABSGULATE TdE UNPERTERBED ORBIT ORBITStO 
0581452 


a: nus SUSSODTINE CALLS Tik FOLLOWING SUBROUTINES: 08811530 
* MORELI = CALCULATE. NEW ELEMENTS AFTER TIME STEP 05811840 
p NEWPOS CALCULATE NEW POSITION AFTER TIME STEP 05811850 
е NE WVEL = CALCULATE NEW VELOCITY AFTER TIME STEE 05811850 
и BIURL 910555 POSIT ICN Омь 08511570 
05311880 

БОЕВЕ PRECISION T, DT, RER, JEP U RIRI IRE RIVI MENR, Y, 05511890 

F PAM, PI A, EOS, TA, Pee ЋЕ БА ЕЕ, ТЕ 0кВ11900 
GRB11910 

0), 03311920 
506) 05811950 
02811940 


(1) 


DIMENSION 22006 UE 
= hee (500) А 


т ШЕШ ТҮСЕ ANOMALY 10 NEGAT Pig SO'ZEAOOP2 «САМЕТ Е ЕСЕ 0) С5811950 
NN US Ci. 6. 21) 5 05811960 

Пан = TA - (22:21) СЕВ11970 

ENDIF СХЕ1 1960 

ORB11990 

* ВЕ 5508р QKkETI; FILL CLOSS TO PERICE 0К812000 

ШКОК (СОСТА „де. 6.221) .&WND. CP. RE: PER)) THEN 05512010 
Oxb120£8 

7: Increment IRUE TINE 05812030 

о БЕ @RB120-0 

КЕМ М LTC AAE, EA, TA, T=, Di, PI, PER) ОКЕ 12050 


CALL NPOS(R1I,RJ,RK,R,LAN,AP,I, TA, ASE 08812060 

CALL NVEL(E,P,TA,LAM, AP,I,VI,VJ, VR, V, MU) 08812070 

08812080 

* INCREMENT STEP COUNTER AND STCRE VALUES CRB12090 
NUM = NUM + 1 ORB12100 

GAEL STORECRI,RJ;RE;R,TA,RIRAY,RJRAÁYjRKRAY, CRB12110 

+ RARAY , TARAY ,NUM, 1, AP, AINRAY, APRAY, ORB12120 

E TT, TIMRAY) ORB12130 
ORB12140 


~ INGREME ND TIME SIEP COUNTER ОКВ12150 

ek + DI ОКЬ12160 

GOTO 230 ORBI2170 

ENDIF 05812180 

RETURN ОКЕ К? ТОО 

END OK512200 

0R512210 

veteveze eese eee eee eere СГ С СК СК ТИТУЛ ORB12220 

ZALCULATE NICSUNDERTURBED .NEWSELEXENTS GRE 12250 
ОО УНА л ПОТЕЗУ УУ i a e a 5с705017 ОХЕ 122 /. 

CRe 2250 


SUBROUTINE NEWELTCMM МА БЕА, ТЕ. ТЕ. В ВЕЕ 
о THIS SUBROUTINE CALCULATES THE Unperturbed NEW ELEMENTS 


THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES: 
NEA = NEW ECCENTRIC ANOMALY 
NTA = NEW TRUE ANOMALY 


DOUBLE PRECISION MM MBSE RES. A, ЕРЕН 


* Increment TIME OF FLIGHT AND CHECK IF TF GREATER THAN PERIOD 
TF = TF + DT 
IF (TF .GT. PER) THEN 
TF = TF - PER 
ENDIF 


* CALCULATE MEAN ANOMALY AND USE TO FIND ECCENTRIC Anomaly THEN NEW 
х TRUE ANOMALY 

MA = MM(TF) 

CALL NEA(MA,E,EA) 

CALL NTA(EA,E,TA,PI) 

RETURN 

END 


jesesevesezevesevedesestevevesevedeseste ede JeJeses'ee de e esee ve еж esee evevesfe ve ese ee ey ede esee e У Је је је ге 


CALCULATE PERTURBED ORBIT 


sesesesevesevevesedeseveves eese eye у о еее еее еее еее еее ее еее еее еее еее еее 


SUBROUTINE -PRETUR( DT, PER AL LAN AR T RERI KENE 
+ - УГ. УЗ У V ВК Ез FW MO PISH ASE в Е 
+ TF,T,NUM,RIRAY,RJRAY,RKRAY, RARAY, TARAY, AINRAY, APRAY,TIMRAY, 
+ TYT,TFEA,TFSU,TFMO,TFDRA,TDI,TDA, IDE, TDMM, TOMA, TDLAN, [DHS Е 


THIS SUBROUTINE CALCULATES THE=SPERUCREED КЕЛЕ 

THIS SUBROUTINE CALLS THE FOLLOWING SUBROUTINES: 

TFORCE = CALCULATE THE TOTAL PERTURBING FORCE ON THE SATELLITE 
PNEWEL = CALCULATE THE Perturbed NEW ELEMENTS 

NPOS = NEW POSITION AFTER TIME STEP 

NVEL = КЕЙ УЕТОСТТҮ АЕТЕК ТТЕ ЛЕР 

PERIOD = PERIOD OF PERTURBEDTORETI 

STORE = STORE POSITION AND ELEMENTS IN ARRAYS FOR PLOTTING 


DOUBLE PRECISION T,DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK,V, 
+ FR,FS,FW,MU,PI,H,A,E,N,TA,P,MM,MA,EA,TF,TT, 

+ — DI,DA,DE,DMM,DMA,DLAN,DH,DAP,EI,EJ,EK, HI, HJ, LP, M, 

+ — DVR,DVS,DVW,DVI,DVJ,DVK 


DIMENSION  RARAY(500),TARAY(500),RIRAY(500),RJRAY(500), 
T RKRAY( 500) , AINRAY( 500) , APRAY( 500) , TIMRAY( 500) 


* SET MEAN RADIUS OF EARTH 
RE = 6400.0 


DT = PER/50 

T = DT 

LE ТА. брате ЈУ ales 
TA = TA - (2*PI) 
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ORB12260 
ORB12270 
ORB12280 
ORB12290 
ORB12300 
ORB12310 
ORB12320 
ORBaZ 330 
ORB12340 
ОКВ12550 
ORB12360 
ORB12370 
ORB12380 
ORB 12390 
ORB12400 
ORB12410 
ORB12420 
ORB 12430 
ORB12440 
ORB12450 
ORB12460 
ORB12470 
ORB12480 
ORB12490 
ORB12500 
ORB12510 
ORB12520 
ORB12530 
ORB12540 
ORB12550 
ORB12560 
ORB12570 
ORB12580 
ОКЕ 590 
ОКВ12600 
0КВ12610 
ORB12620 
ORB12630 
ORB12640 
ORB12650 
ORB12660 
ORB12670 
ORB12680 
ORB12690 
ORB12700 
ORB12710 
ORB12720 
ORB12730 
ORB12740 
ORB12750 
ORB12760 
ORB 12770 
ORB12780 
ОКЕ 12790 
ORB12800 
ORB12810 


+ + + + 


+ + 


241 


ENDIF 
ieee ОЕ PER) THEN 


ИРИНЕ 


EDI 


CONTINUE Around ORBIT FOR ONE PERIOD 
Pee irs. Ul. РЕК) АМР. 


(А РА АРОН НОУ 


INCREMENT TRUE TIME 
TT = TT + DT 
CALL TFORCE(AL,LAN,AP,1,RI,RJ,RK,R,V1,VJ,VK,V, 
TT,FR,FS,FW,MU,PI, 
FEA,FSU,FMO,FDRA,FOR, 
EI,EJ,EK,E,A,T,LP, TA, PER, EA, MA, TF,P, 
MM,N,H,HI,HJ,DT) 
CALL PNEWEL(FR,FS,FW,H,R,A,E,N,TA,DT,1,LAN,AL, 
AP,P,MM,MA,EA,TF,T,MU,PI, 
DI,DA,DE,DMM,DMA,DLAN,DH, DAP) 
CALL NPOS(RI,RJ,RK,R,LAN,AP,I, TA,A,E) 
CALL NVEL(E,P,TA,LAN,AP,1I,V1,VJ, VK, V,MU) 


CALCULATE NEW PERIOD AND RESET TIME STEP AND TIME COUNTER 


IF NOT AT END OF ORBIT 

IF (T .LT. (PER-DT)) THEN 
CALL PERIOD(A,PER,PI,MU) 
ОТ = РЕВ/50 
E 

ENDIF 


INCREMENT STEP COUNTER 
NUM = NUM + 1 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY, 


RARAY,TARAY,NUM,I,AP,AINRAY,APRAY, 


TL, FIRAT) 


TOTAL ELEMENT CHANGES 


TDI = TDI + SNGLCABS(DI)) 
TDA = TDA * SNGL(ABS(DA)) 
ІШЕ = TIDE + SNGLCABS(DE) ) 


ЮВ = TDM + SNGLCABS( DMM) ) 
TDMA = TDMA + SNGLCABS(CDMA) ) 
TDLAN = TDLAN * SNGL(ABS(DLAN)) 
ТОН = ТОН + SNGL(ABS(DH)) 


TDAP = TDAP + SNGLCABS( DAP) ) 
ТЕЕА = TFEA + FEA 
IESU = IFS5U f ESU 
TEO = ТЕМО + FMO 


TFDRA = TFDRA + FDRA 


CHECK FOR IMPACT 

DER ШЕ Гы) ТЕ 
БОЛАП SATELLITE WILL INPACT THE EARTH! 
Т < РЕК 

ENDIE 


TC RESES TATIE COUNTER 
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ORB12820 
ORB12830 
CRB12840 
ОКВ1 2850 
ОКВ12860 
ORB12870 
ORB12880 
06812530 
ОКВ1 2900 
ORB12910 
ШЕВІ2920 
0КВ12930 
ORB12940 
ORB12950 
ORB12960 
98812520 
ОКВ12980 
КЕ Ша 20) 
ОКВ13000 
ORB13010 
ORB13020 
ORB IS030 
ORB13040 
ORB13050 
ORB13060 
ORB13070 
ORB13080 
ORB13090 
ORB IS100 
ORB13110 
ORB13120 
ОВ) 
ОКВ13140 
ORB13150 
ORBIS T60 
ORB13170 
ORB13180 
ORB13190 
ORB13200 
ORB13210 
ORB13220 
ORB13230 
ORB13240 
ORB13250 
ORB13260 
Во 
ORB13280 
ORBI 3290 
ORB13300 
ORB13310 
(ВО О 
ORB 13350 
ORB13340 
ORB13350 
ORB13360 
ORB13370 


= Ре О 
GOTO 240 
вв 
RETURN 
END 


yes eyesevev ese eee ese des eee eve eee ev evesfedese eee es eese ees eese eye ye ee vec eese eese eve sees eee vestes 


* CALCULATE THE PERTURBING FORCES 
Jes ees eyes eese ev eye А И esee ey eee 


SUBROUTINE ТЕОКСЕСАЬ, ВАМ, АР, КОК ТЕ 


ar ER,ES,EW,MU,PISBEASESUSBISMESRASFOR, 

* EI EJ EK,E,A T LP, ITAS FERSEN MA IF P, 

a MNA HIKE DID 
THIS SUBROUTINE SUMS ALL THE PERTURBING FORCES FOR THE TOTAL 
P PERDUBSBPNGODORCES 


THE FOLLOWING SUBROUTINES WERE CALLED: 


OBERT = OBLATENESS OF THE EARTH 

FSUN = GRAVITATIONAL Attraction OF THE SUN 
* FMOON = GRAVITATIONAL Attraction OF THE MOON 
э FDRAG = DRAG FORCES 


DOUBLE PRECISION FER,FES,FEW,FSR,FSS,FSW,FMR,FMS,FMW,MU,PI, 
+  FDR,FDS,FDW,FR,FS,FW,RI,RJ,RK,R,AL,1I,TT,LAN,AP,VI,VJ,VK,V, 
+ EI,EJ,EK,E,A, T, LP, TA, PER, EA, MÀ, TF,P, 

E MM,N,H,HI,HJ,DT 


CALL OBEART(RI,RJ,RK,R,AL,I,FER,FES,FEW,MU) 

CALL FSUN(TT,RI,RJ,RK,R,FSR,FSS,FSW,PI) 

CALL FMOON(TT,RI,RJ,RK,R,FMR,FMS,FMW,PI) 

CALL FDRAG(RI,RJ,RK,R,VI,VJ,VK,V, LAN, AP, I, FDR,FDS,FDW, 

+ EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 
+ MM,N,H,HI,HJ,DT) 


SUM VECTOR FORCES 


ЕК = ЕЕК + ЕТ 5 ЕЕЕ 
FS = FES + Е ВИ ЕО - 
Fw = FEW + FSW + FMW + FDW 


CALCULATE TOTAL FORCE FROM EACH, AND TOTAL OF ALL 


FEA = SNGLCSQRT( (FER**2)+( FES***2 0) ЕЕС 2000) 
FSU = SNGLCSQRT( (FSR**2)+(FSS**2)+( FSW**2) )) 
ГМО = SNGLCSQRT( (FMR**2)+(C FMS**2)+( FMW**2) ) ) 


FDRA = SNGL(SQRT( (FDR***2 )+( FDS****2 )+( FDW?*2 ) ) ) 
FOR = SNGLCSQRT( (CFR**2)+(C FS**2)+( FW**2) ) ) 


RETURN 
END 


оу јео у УУ Је У де С МУ у о С КСО ПСК СЗ ОКС ССЗ о У пио 
SUBROUTINE ОВЕАСТ ЕКТ КИ КК KR An PER. FES eee ion 


= THISGSSUBROUTINE GADCUBATES THESPEERTURDTIAGOEIDSISCEMEDUESSPOTME 
т OPLIQVENESS OF THE BARTE. 


ORB13380 
ORB T3320 
ORB13400 
ORB13410 
ORB13420 
ORB13430 
ORB13440 
ORB13450 
ОКВ13460 
ОКВ13А 70 
ORB13480 
ORB13490 
ORB13500 
ORB13510 
ORB13520 
ORB13530 
ORB13540 
ORB13550 
ORB13560 
ORB13570 
ORB13580 
ORB13590 
ORB13600 
ORB13610 
ORB13620 
ORB13630 
ORB13640 
ORB13650 
ORB13660 
ORB13670 
ORB13680 
ORB13690 
ORB13700 
ORB13710 
ORB13720 
ORB13730 
ORB13740 
ORB13750 
ORB13760 
ORB13770 
ORB13780 
ORB13790 
ORB13800 
ORB13810 
ORB13820 
ORB13830 
ORB13840 
ORB13850 
ORB13860 
ORB13870 
ORB13880 
ORB13890 
ORB13900 
ОКВ1ЗО10 
ОКВ13920 
ОКВ13930 


ВРЕТ ПРОВЕО 272 RE FER FES FEW RI; RJ RRR AL ILMU 


72 = 1.082364D-03 
mee = б 37620703 
ПИ (( -3. 0DT00*MU*J2*(RE-*225)7(2- OD*00*CR**49)955* 
F OD OORE ОЗ ОООО СООО S2) * CUODSTINCAT))7*2))) 
FES = ((-3. 0OD*00* MU* J2* ( RE**2))/ (R***4))* 
+ (CCDSINCI))**2)* (DSINCAL))*(DCOS(AL))) 
FEW 2 ((-3. 0D*00*MU* J2* ( RE***2) ) / (P054) )** 
t (DSINCI)*DCOS(I)*DSIN(AL)) 
RETURN 
END 


sce see ede ses Aes es eseves eee eve УУ их УУ a 4 


SEESESOUTINECESUNCITRI,RJ,RK,R;FSR,FSS,FSW,PI) 


* Mims SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO THE SUN 
* THE FOLLOWING SUBROUTINES ARE CALLED: 
ve SUNPOS SUNS POSITION ORBITING AROUND EARTH 


7% HEVBOD PERTURBING FORCE FROM A Heavenly BODY 
POELE PRECISION PSReFSS¢ row, 1. 
Es LEM ESJSSK SDEANSSISSABSSHUSTT,RIJRJJRK,R,RS 


* SUNS GRAVITATIONAL PARAMETER 
Bae 1, 32715440+11 
CALL SUNPOS(TT,RSI,RSJ,RSK,RS,SLAN,SI,SAL,PI) 


ORB13940 
ORB13950 
ORB13960 
ORB13970 
ORB13980 
ORB Tosa” 
ORB14000 
ORB14010 
ORB14020 
ORB14030 
ORB14040 
ORB14050 
ORB14060 
ORB14070 
ORB14080 
ORB14090 
ORB14100 
ORB14110 
ORB14120 
ORB14130 
ORB14140 
ORB14150 
ORB14160 
ORB14170 
ORB14180 
ORB14190 
ORB14200 
ORB14210 
ORB14220 


CALL HEVBOD(RI,RJ,RK,R,RSI,RSJ,RSK,RS,SLAN, SAL, SI ,SMU,FSR,FSS,FSW)ORB14230 


RETURN 
END 


о о о о ыы ыы об anis ef auto 
SR 20 Ce Oe 72V UV VV VV 2V V8 VV EN GV EN OV VV UY ON V» 2& QV OV VV VV VV Ov 2V OV GN ON ИУ 2% 2% 2% 2% 2% 2% 24 2% UV EV V* 2v v VV VC VV 2V UN (зг 


EESEDUTINE FHOONCTT,RI,RJ,RK,R,FMR,FMS,EMW,PI) 
THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO The MOON 


* THE FOLLOWING SUBROUTINE ARE CALLED: 
у MONPOS = MOONS POSITION ORBITING AROUND THE EARTH 
т: rev BOD = PERTURBING FORCE FROM A HEAVENLY BODY 


ROUBLE PRECISION FMR,FMS,FMW,RMI,RMJ,RMK,MLAN,MI,MAL,MMU, 
J TERT RIRES RRM EI 


х MOONS GRAVITATIONAL PARAMETER 
MMU = 4. 90287D+03 


ша НОМРОЗСТТ, ВИТ, ЕМУ, КМК ВЫ, МЬАМ  М:,МАБ,РТ) 


ORB14240 
ORB14250 
ORB14260 
ORB14270 
ORB14280 
ORB14290 
ORB14300 
ORB14310 
ORB14320 
ORB14330 
ORB14340 
ORB14550 
ORB14360 
ORB14370 
ORB14380 
ORB14390 
ORB14400 
ORB14410 
ORB14420 


МИР ИО ВОДИО БТ КК Е КР КУУ КНК, КМ, МГАМ МАТ, МТ ММ ЕМК,„ЕМ5, ЕМИ) ОКВТААЗО 


RETURN 
END 
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ОКВ14440 
ORB14450 
ORB14460 
ORB14470 
ORB14480 
ORB14490 


sc 
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ela cla nl a ъв? въ? йә ъ! э ә ы? а ъ! ә ьа ъЇ ә ъ а ъ' э ъ! ә ъ а ўа ъ! э ә ъ э ы? ә ы! чы! ә ъ ! э чы] ә чы! ә ъЇ э ъй тъ? а «8, У У ufa nf a nf nf afa eta e у: > ul o nf ~! јаја ја одаја ~! ufo afa nda afa n anf nf anf 
0% 02% ғ% FR FR PR Pd Fd ER. ER EL ER FR FR 04V OV 68 4 0S 0S 04S 9S9 съ 9S 9S 0S 06S съ 4ъ 0:09 2% 9% 4% vee. еъ оъ тъ еъ venr % 2% 2. ou ee 4% 727% съ 4% % 4. 22%: 0% 9% 9% Фу 27% e^ 


a's 
* 


ж 


T 
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ВИК ЕНӘ ЕН) 
THIS SUBROUTINE CALCULATES THE PERTURBING FORCE DUE TO A 
HEAVENLY BODY. 


THE FOLLOWING SUBROUTINE WAS CALLED: 
IJKRSW = 'IJK' SYSTEM TO THE 'RSW' SYSTEM 


DOUBLE PRECISION DOTS PRIS FHI EAE REE ГОРЕ ere ep Ер 


LAN,AL,INC,MUP,I,J,K,IP,JP,KP,M,FHR,FHS,FHW 


ORB14500 
ORB14510 
ORB14520 
ORB14530 
ORB14540 
ORB14550 
ORB14560 
ORB14570 
ORB14580 
ORB14590 


CALCULATE UNIT VECTOR FOR SATELLITE AND PERTURBING BODIES POSITIONORB14600 


I = RI/R 
RITR 
RK/R 
= eRe ly RE 
JP = RPJ/RP 
= RPK/RP 


Са 
ІШІ 


CALCULATE DOT PRODUCT OF UNIT VECTORS 
DOT = (( I*IP 1 ТОТ 


CALCULATE FORCES IN THE 'IJK' SYSTEM 


ЕНІ = (MUP/(RP**2))*(R/RP)*(3. 0D*00*DOT*( IP) -(I)) 
ЕНТ = (MUP/(RP**2))-(R/RE) (32 0Dt00DOl 3 a) 
FHK = (MUP/(CRP**2))*(R/RP)*(C 3, 0D+00- DOr ei) 


Transform FORCES TO THE RSW SYSTEM 

CALL IJKRSW(LAN,AL,INC,FHI,FHJ,FHK,FHR,FHS,FHW) 
RETURN 

END 


SUBROUTINE SUNPOS(TT,RSI,RSJ,RSK,RS,SLAN,SI,SAL,PI) 
THIS SUBROUTINE CALCULATES THE SUNS POSITION 


VARIABLES USED TO DESCRIBE THE SUNS ORBIT: 
SI = SUNS INCLINATION 
SLAN= SUNS Longitude OF ASCENDING NODE 


SAP = SUNS ARGUMENT OF PERIGEE 
RS = SUNS ORBITAL RADIUS 

STA = SUNS TRUE ANOMALY 

SAL = SUNS ARGUMENT OF LONGITUDE 


DOUBLE PRECISION SLAN,SI,SAL,RS,STA,SAP,TT,RSI,RSK, 


+  RSJ,RSP,RSQ,RSW,PI 


SI = 4.09279709D-01 

SLAN = 0. 0р+00 

ЗАР = 0. 0р+00 

RS = 1. 4959965D+08 

STA = ((2. O*PI)/( 365.0 * 86400.0) * TT) 
SAL = STA + SAP 


CALCULATE SUNS POSITION IN ‘PQW' SYSTEM 
RSP = RS*DCOS(STA) - 


ORB14610 
ORB14620 
ORB14630 
ORB14640 
ORB14650 
ORB14660 
ORB14670 
ORB14680 
ORB14690 
ORB14700 
ORB14710 
ORB14720 
ORB14730 
ORB14740 
ORB14750 
ORB14760 
ORB14770 
ORB14780 
ORB14790 
ORB14800 
ORB14810 
ORB14820 
ORB14830 
ORB14840 
ORB14850 
ORB14860 
ORB14870 
ORB14880 
ORB14890 
ORB14900 
ORB14910 
ORB14920 
ORB14930 
ORB14940 
ORB14950 
ОКВ14960 
ОКВ14970 
ОКВ14980 
ОКВ14990 

ОКВ 15000 
ORB15010 
ОКВ15020 
ОКВ15030 
ORB15040 
ORB15050 


RSQ 
RSW 


RS*DSIN(STA) 
0. 0D«00 


ШІ 


TRANSFORM POSITION TO 'IJK' SYSTEM 

CALL PQWIJKR(SLAN,SAP,SI,RSP,RSQ,RSW,RSI,RSJ,RSK) 
RETURN 

END 
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Е 


ео)" 
0% 9% 60% дъ 6*9 06^ 


BEENOUTPNESMONPOSCET,RAI,RAJ,REHR, RM, MLDAN, MI , MAL, PI) 
THIS SUBROUTINE CALCULATES THE MOONS POSITION 


ШЕТАБЕЕ5 USED TO DESCRIBE THE SUNS ORBIT: 
MI = MOONS INCLINATION 
MLAN= MOONS Longitude OF ASCENDING NODE 


MAP = MOONS ARGUMENT OF PERIGEE 
RM = MOONS ORBITAL RADIUS 

MTA = MOONS TRUE ANOMALY 

MAL = MOONS ARGUMENT OF LONGITUDE 


WOU BER PRECISION MI,MLAN,MAL,RM,IM,MTA, RMP ,RMO,RMW , 
ГИБКО МАРЗ РІ 


МІ = 4.991651660-01 
RM = 3.844D+05 


MLAN = 0.0 

ENES- 0002. 0-PI)/(27.3 * 3600) ^ TT) 
МАР = 0. 00+00 

MAL = МТА 


CALCULATE MOON POSITION IN ‘PQW' SYSTEM 
КІР RM**DCOS( MTA) 

RQ RM*DSIN(MTA) 

КМУ 0 


TRANSFORM POSITION TO ‘IJK’ SYSTEM 

ERER POK TIKLAN, MAP MI, RMP, RMO, RW, RMI, RMJ, RMK) 
RETURN 

END 


BRON aN VV ev 4 dv VEN eV 2v 4S EV YN Ov VV EN VV Qv aV VV UV Ev VV QV EV VV ON QN EV eV VN ON UN OV Ov VV ev EV Ov av ev VA Ev ev 0 Uv av 2v Vk Vv 2v ON Av VR EN EN VV VV VN 


SUBRROUTTNESCEBRAGCRI.RJ,RK,R,VISVJS;VR,V, DAN, AP, T; FDR, FDS,FDVW, 


+ EI,EJ,EK,E,A,T,LP,TA,PER,EA,MA,AL,TF,P,PI,MU, 
т MM,N,H,HI,HJ,DT) 


mils SUBROUDINE  GALGULATES THE PERTURBING FORCE DUE TO DRAG 


TEE FOLLOWING VARIABLES ARE USED TO MODEL THE ATMOSPHERE: 


RE = RADIUS OF EARTH 

M = MASS OF SATELLITE 

AR = FRONTAL SURFACE AREA OF SATELLITE 
Z =ОШ ОЕ OF SATELLITE 

К = EXPONENTIAL DECAY FACTOR 

BESO = ORAL DENS. Ty 

C2 UOSERISIEVIESUT DES 


fA 
„ә 


ОКВ15060 
ORB15070 
ORB15080 
ORB15090 
ORB15100 
ORB15110 
ORB15120 
ORB15130 
ORB15140 
ORB15150 
ORB15160 
ORB15170 
ORB15180 
ORB15190 
ORB15200 
ORB15210 
ORB15220 
ORB15230 
ORB15240 
ORBI5250 
ORB15260 
ORB15270 
ORB15280 
ORB15290 
ORB15300 
ORB15310 
ORB15320 
ORBIS 330 
ORB15340 
ORB15350 
ORB15360 
ORB15370 
ORB15380 
ORB 5290 
ORB15400 
ORB15410 
ORB15420 
ORB15430 
ORB15440 
ORB15450 
ORB15460 
ORB15470 
ORB15480 
ORB15490 
ОКВ15500 
ORB SS TIO 
0К815520 
ORE 135.30 
ORB15540 
ORBISSSD 
ORB15560 
ORB15570 
ORB15580 
ORB15590 
ORB15000 
ОКВ15610 


Ааа ОН УУ УУ И У: 


+ CALCULATE PERTURBED ud HORE 


DOUBLE PRECISION MAG,M,K,FDR,FDS,FDW,RE,AR,Z,DENO,CD,DEN, 
: FDJ,FDR,FDI,RI,RJ,RK.VI.VUS VR V. DAN АР, PARS 
EI,EJ.EX,E. A, T, LP, TA, PER, EA, MÀ, AL, TF, P, PI,MC, 


Мә ПО ЈА А У КО У РАИ ОР Јн 


КЕ - 622781830712 


| = 1.00702 
АК = 2. 010-01 
д. = В ВЕ 


DEPENDING ON ALTITUDE SET ATMOSPHERE VARIABLES 
ТЕСЕ БОЕО Е 

К = 4. 740-02 

DENDO 1.225050 

Ср = 1. 00+00 
ЕБЕТ. СЕ 5 УПАО JS IHES 

К = 3.46140-02 


DENO — 17995569209] 
С) = 2. 08-95 
ЕВЕ 


К = 2. 21636056 
= 1.015484D-07 
CD = 2. 90595 


CALCULATE ATMOSPHERIC DENSITY 
DEN >= DENO ВЕХРОБЕ А) 


CALCULATE MAGNITUDE OF DRAG FORGE CAND CITAT РОА ИЕ и 
MAG = -(0. 5D+00)=CD=AR= DEN Ya ( 1. 0D-05)/M 
ЈЕ (АВУСМАВ) „БТ. #20 АТИКА 
Мас - -1- 90-29 
ENDID 


GIVE DRAG FORCE A Direction Orm IAC SS ER EEOC 


Ррк- 0.0 

FDS = MAG * V 
Dw - 00 

RETURA 

END 


Јеху по 
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SUBROUTINE PNEWEL(FR,FS,FW,H,R,A,E,N, TA,DT, I, LAN, AL, AP, P, 

+ MM,MA,EA.TE, T, MU, PI,DI,DA, DE, DMM,DMA , DLAN, DH,DAP) 

THIS SUBROUTINE CALCULATES THE NEW ELEMENTS FROM THE PREVIOUS 
ELEMENTS ADDED TO THE RATES OF CHANGE FOR ONE STEP 


THE FOLLOWING SUBROUTINES ARE CALLED: 


RATE = CALCULATES RATES OF CHANGE OF ORBITAL ELEMENTS 
NANGMO = NEW ANGULAR MOMENT US DERN) 
мола = NEW SEVIS ТОК AXIS (Ори) 


wrap 
ЧАО и 


мм, ВСЕ TRICITY ( NEWE) 


ORB15620 
ORB15630 
ORB15640 
ORB15650 
ORB15660 
ORB15670 
ORB15680 
ORB15690 
ORB15700 
ORB15710 
ORB15720 
08815736 
ОКВ15 740 
ORB15750 
ORB15760 
ORB15770 
ORB15780 
ORB15790 
ORB15800 
ORB15810 
ORB15820 
ORB15830 
ORB15840 
ORB15850 
ORB15860 
ORB15870 
ORB15880 
ORB15890 
ORB15900 
ORB15910 
ORB15920 
Око о 
ORB15940 
ORB15950 
ORB15960 
ORB15970 
ORB15980 
ORBIS 72 ¢ 
ORB16000 
ORB16010 
ORB16020 
ORB16030 
ORB16040 
ORB16050 
ORB16060 
ORB16070 
ORB16080 
ORB16090 
ORB16100 
ORB16110 
ORB16120 
ORB16130 
ORB16140 
ORB16150 
ORB16160 
ORB16170 


( NEWAP) 


Е НОЩА ТОТ ве АТ AP ВО 


NINCL = NEW INCLINATION CNEWT) 
NASNOD = NEW LONGITUDE OF ASCENDING NODE (NEWLAN) 
Кс ГЕК = New ARGUMENT OF PERIGEE 
NMNMO = у MEAN MOTION (КЕМ) 

О = ова ОГО CIL) 

NMMAN = NEW MEAN ANONALY (NEWMA) 
NEA = NEW ECCENTRIC ANOMALY (EA) 
NTA = NEW TRUE ANOMALY (TA) 
мие = TIME OF FLIGHT (TF) 

ШЕРГЕ PRECISION  FR,FS,FW,DMM,H 
и ММ 


+ NEWH,NEWA,NEWE 


SIMEACLTP T. MUSPISDASDHSJDESDISDDANSDSP.DMA: 


, NEWI , NEWLAN , NEWAP , NEWMM 


INCREMENT TIME OF FLIGHT BY ONE TIME STEP AND CALCULATE RATES 


TF = 


ШЕБЕРІ 


ШЕШІ KATBESCDH,DA,DE,.DI,DLAN;,DAP,DMM,DMA,E,; 1, R,A,FR,ES,FW, 


~ 


ТАСАГ Н, РИКИ) 


CALCULATE NEW ELEMENTS 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


SET 
A 


L 


Í 


нии 


ПЕСО ООН, МЕНН) 
NSMACA ,DT,DA,NEWA) 

NMEGCUE DT DESNEWEJ 
ДОС РА ЉЕ ХЕМЕ) 
NASNOD( LAN, DT, DLAN , NEWLAN ) 
NARPER( AP ,DT, DAP, NEWAP) 


ЕШ ШЕ Т5 TO NEW ELEMENTS 
NEW aA 
NENE 


RESI 


LAN = NEWLAN 


AP = 


РЕДА: 


КОЕ 
CALL 
CALL 
CALL 
CAL 
CALL 
AL = 


DUE 
и шев ар) 


ШОО ТЕО ТЕ ONS TIME STER 
ПА МО АМА НО) 

ПА МА ВОТ, ТЕ UNA, PT) 
NEACHA,E,EA) 

КОЛГЕ ЕТА УРТ.) 
ШЕШН ма ТЕ) 

ТА + АР 


RETURN 


END 
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SUBROUTINE KRATESCDH .DA,DE,DI,DLAN, DAP, DUM, ,DMA,E,MM,R,A,FR,FS,FW, 


ey Us 


Ew DERE T) 


ЕН ШӘ SUBROUTINE Calls THE FOLLOWING SUBROUTINES TO CALCULATE THE 
Di вв СЕ OF THE ORBITAL ELEMENTS: 


Эа 


Та а 
КЕС 
г 
AAA a NN 


RATE-Or -CHANGE 
BATESOE-GDEANGE 
AATZ-UOr-CHANGE 


OF 
OF 
Ог 


TRE 
THE 


53 


THE SE 


ПЕ ВА 15 14) 


ECCENTRICITY KDE) 
INCLINATION (DI) 


ORB16180 
ORB16190 
ORB16200 
ORB16210 
ORB16220 
ОБЕ ОР 30 
ОКВ16240 
ОКВ16250 
ОКВ16260 
ORB16270 
ORB16280 
ОВ ево 
ОКВ16300 
0КВ16310 
ORB16320 
ORB T6320 
ОКВ 163420 
ОКВ16350 
ОКВ16360 
ORB16370 
ORB16380 
ORB16390 
ORB16400 
ORB16410 
ORB16420 
ORB16430 
ORB16440 
ORB16450 
ORB16460 
ORB16470 
ORB16480 
ORB16490 
ORB16500 
ORB16510 
ORB16520 
ORB16530 
ORB16540 
ORB16550 
ORB16560 
ORB16570 
ORB16580 
ORB16590 
ORB16600 
ORB16610 
ORB16620 
ОКВ16630 
ОКВ16640 
ОКВ16650 
ОКВ16660 
ОКВ16670 
ORB16680 
ORB16690 
ORB16700 
ОКВ16710 
ОКВ16720 
ORB10730 


+ 


' ? 
песен 


RLAN = RATE-OF-CHANGE OF THE Longitude OF THE ASCENDING NODE 
(DLAN) 

RAP = RATE-OF-CHANGE OF THE ARGUMENT ОЕ БЕРЕР СТАР) 

nol = RATE-OF-Change OF THE MEAN ЧОТТОУ (0) 

RMA = RATE-OF-CHANGE OF THE MEAN ANOMALY (DMA) 

RANGMO = RATE-OF-CHANGE OF THE ANGULAR MOMENTUM (DH) 

DOUBLE PRECISION DH,DA,DE,DI,DLAN,DAP | DMM ОМА ЕМИ КО АЕК БОНИ 
TAAAC H Е В HP 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


RSMAX(E,MM,R,A,FR,FS,DA,TA) 
RECC(E,MM,R,A,FR,FS,TA,DE) 
RINC(E,MM,R,A,FW,AL,DI) 
RLAN(E, MM, R, A, I,FW, AL, DLAN) 
RAP(E,MM,R,A,I,H,P,AL,TA,FR,FS,FW,DAP) 
RMM(MM,A,DMM,DA,MU) 

CALL RMA(E,MM,R, A, TA, DMM,FR,FS,DMA,T) 

CALL RANGMO(R,FS,FW,DH) 

RETURN 

END 


о 
ЗЕВКОШТТКЕ ВАО В Е.В ВЫ) 

THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE 
ANGULAR MOMENTUM 


DOUBLE. PRECISION MES FPW DEW DESS ONTE 


DHW = К * FS 

DHS = R * Fw 

DH = DSQRTCCDHW***2) + CDHS***2) ) 
RETURN 

END 
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SUBROUBINE RSMAXCE ВОК. А PR TS DANTA) 
THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE SENT ASTR 
AXIS 


DOUBLE PRECISION DA,FR,FS,E,MM,R,A,TA,ET 


TRAP (E) SO DENOMINATOR DOES NOT GOTO ZERO 
IF (E -GT OTO ATEEN 


El gm 
ELSE 

ЕТ + Е 
ENDIF 


DA = ((2. OD+00"E “DSINCTA D PI = DSORIC TOPOR E ОНИЕ 


+  ((2. 0D+00*A*DSQRT(1. 0D+00-(E **2)))/(MM *R))*FS 


+ . . 
ААА ААА 


RETURN 
2556, 


.: 52-4 ЖЕ пере но в ров го про в па MES NONE S MI RE 
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ORB16740 
ORB16750 
ORB16760 
ORB16770 
ORB16780 
ORB16790 
ORB16800 
ORB16810 
ORB16820 
ORB16830 
ORB16840 
ORB16850 
ORB16860 
ORB16870 
ORB16880 
ORB16890 
ORB16900 
ORB16910 
ORB16920 
ORB16930 
ORB16940 
ORB16950 
ORB16960 
ORB16970 
ORB16980 
ORB16990 
ORB17000 
ORB17010 
ORB17020 
ORB17030 
ORB17040 
ORB17050 
ORB17060 
ORB17070 
ORB17080 
ORB17090 
ORB17100 
ORB17110 
ORB17120 
94:27) 11570. 
ORB17140 
ORB17150 
ORB17160 
ORB17170 
ORB17180 
ORB17190 
ORB17200 
ORB17210 
ORBIT2260 
ORB17230 
ORB17240 
ORB17250 
ORB17260 
ORB17270 
ORB17280 
ORB17290 


а ге Гател 


”, 


+ 


ООСС ТВ КЕСССЕ ММС ВОЈА РЕКЕ Ео ТА БЕО 


THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE ECCENTRICITY 


ЕЕ КЕ Т5ТОХ DESER: FS E МОК. А,ТАЈЕТ 


TRAP (E) SO DENOMINATOR DOES NOT GOTO ZERO 


IF (E.LT. 0. 1) THEN 
ЕТ - 0.1 

ELSE 
ET=E 

ENDIF 

DE = ((DSQRT(1. OD+00 - (E %**2))*SIN(TA))/(MM *A))*FR + 
((DSQRT(1. OD+00 - (Е ##2)))/(ММ *ET*( A%*2)))* 
((А%%2У5(1 00400 - (Е %527)/(Е) - (ЕУ)%Е5 

RETURN 

END 


t 


t 
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SUBROUTINE RLANCE ,MM,R,A,1I,FW,AL, DLAN) 


THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE LONGITUDE 
OF THE ASCENDING NODE 
DOUBLE PRECISION DLAN,FW,E,MM,R,A,I,AL,ET,IT 
TRAP (E) AND (I) SO DENOMINATOR DOES NOT GOTO ZERO 
ШЕ (Е СТ. 0. 9) ТНЕЧ 
О 9 
ELSE 
ET =E 
ENDIF 
lect. LT, 0,01745) THEN 
Ша 0.01745 
ELSE 
т 
в В.Р 
DLAN 2 (R*FW*DSIN(AL))/(MM **(A**2)*DSQRT(1. 0D400 - (ET4*2))* 
DSIN(IT)) 
RETURN 


END 
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ШЕПІСІТКЕ ҚКАР(Е,ММ,К,А,І,Н;Р,АІ,,ТА,ҒЕ,Е5,ҒИ,ПАР) 
THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE ARGUMENT 
ШЕ PERIGEE 


B E KEGISTON 2 DAPRK,VAFS, DAPW ПАР ЕК Е Ви Е МОК ЈН Р,АБ,ТА, 


аи ЕЛ А 


ТКАР (Е) AND (I) SO DENOMINATOR DOES NOT GOTO ZERO 


ПИАТ ООП та У THEN 
IT = 0. 01745 

LSE 
I I 


к 
jag 
|| 


А . 4 


tA 


ORB17300 
ORB17310 
СЕРІ 920 
08817330 
ORB17340 
ORB17350 
ORB17360 
ORB17370 
ORB17380 
ORB17390 
ORB17400 
ORB17410 
ORB17420 
ORB17430 
ORB17440 
ORB17450 
ORB17460 
ORB17470 
ORB17480 
ORB17490 
ORB17500 
ORBIZSIO 
ОКР 20 
О 520 
ОКВ1 7540 
ОКВ1 7550 
ОКВ 1 7560 
ОКВ 1 7570 
ОКВ1 7580 
ORB17590 
ORB17600 
ORB17610 
ORB17620 
ORB17630 
ORB17640 
ORB17650 
ORB17660 
ORB17670 
ORB17680 
ORB17690 
ORB17700 
(ОВ 710 
ORB17720 
ORB17730 
ORB17740 
ORB17750 
ORB17760 
ORB17770 
ORB17780 
ORB17790 
ORB17800 
ORB17810 
ORB17820 
ORB17830 
ORB17840 
ORB17850 


шыт О Ол ТЕМ 


ЕТ = 0.9 
ЕЦЗЕТЕ (Е ПЛ. 0.1) SEHE 
ЕР = О 
ENSE 
ЕТ = Е 
ENDIF 
DAPR = (-DSQRT(1.0+00 - (E **2))*DCOS(TA))/(MM *A*ET) * FR 


DAPS = (P/(ET*H))*(DSIN(TA) )* 
(1. 0D+00 + 1. OD+00/(1. OD+00 + ET*DCOS(TA))) *FS 
DAPW = (-R*(1. OD+00/DTAN( IT) )*DSIN(AL) )/ 
(MM *(A**2)*DSQRT(1. 0D+00 - (ET**2)))*FW 
DAP = DAPR + DAPS + DAPW 
RETURN 
END 


еее еее еее ее kekeke kekeeke 


SUBROUTINE RINCCE,MM,R,A,FW,AL,DI) 
THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE INCLINATION 


DOUBLE PRECISION ПТ, ЕМ Е И КОД АНТЕ 


TRAP (E) SO DENOMINATOR DOES NOT GOTO ZERO 
IF (E GT-02) ATHEN 


= Оз 
IPSE 

ЕТ + Е 
ENDIF 


DI = (R“FW*DCOS(AL))/CMM #(A*#2)*DSORT(1. OD+00 - (ET**+2))) 
RETURN 
END 


9v (VIV IN (y 98 4v VV ENED ER ENON ER GLEN GO CC EN OD EN ENED EN ENED EDC TL ER EN FON OL EL ELEN EN FR IL ECTS ELEN VED EVER ER GL ESOL ES OLED EDV FLED ERED GN BL ERED 0% 


SUBROUTINE КЕБЕ У 
THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE MEAN MOTION 


DOUBLE PRECISION DMM,DA,MM,A,MU 

DMM =((-3. OD+O0*MU)/( 2. OD+O0*MM *CA**4)))* DA 
RETURN 

END 


2729622267:72269:9296265:20925692269226222929292929292579296969259292969292929259279292929699296909259969290963696922 09092 


SUBROUTINE КМА(СЕЈУМ,Е,А, ТА DAIM FR ЕЗ БЛ ИЫ 


THIS SUBROUTINE CALCULATES THE RATE OF CHANGE OF THE MEAN Anomaly 


DOUBLE PRECISION DMAA,DMAB,DMAC ,DMAD,DMM,FR,FS,DMA,E,MM,R,A,TA, 


- E 


TRAP (E) SO DENOMINATOR DOES NOT GOTO ZERO 
ТОО ЕЕ сто. оу ты 

ЕТ + 0.9 
Беър (ТОЛТО. Ту ТШЕ: 
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ОКВ17860 
ОКВ17870 
ОКВ1 7880 
ОКВ1 7890 
ОКВ17900 
ОКВа ЛОШО 
ORB17920 
ORB17930 
ORB17940 
ORB17950 
ORB17960 
ORE T7 970 
ORB17980 
ОКВ1 7990 
ORB18000 
ORB18010 
ORB18020 
ORB18030 
ORB18040 
ORB18050 
ORB18060 
ORB18070 
ORB18080 
ORB18090 
OKB18100 
ORB18110 
ORB18120 
ORB18130 
ORB18140 
ORB18150 
ORB18160 
ORB18170 
ORB18180 
ORB18190 
ORB18200 
ORB18210 
ORB18220 
ORB18230 
ORB18240 
ORB18250 
ORB18260 
ORB18270 
ORB18280 
ORB18290 
ORB18300 
ORB18310 
ORB18320 
ORB18330 
ORB18340 
ORB18350 
ORB18360 
ORB18370 
ORB18380 
ORB18390 
ORB18400 
ORB18410 


PT >‘ ~! ~! ој аһ е PU PUPPES PRI ate ate Ses e$. 
еъ съ гъ оу ЕУ гу v9 Co vv 4^5 27% 4% 4% 2% дъ 4% 


- 2 m .7Ҙ... ae "el ont n ae ae о ~! ~ ae m “2... P m m" ~. “..... m m m -. “".-!. alano als m a's пије ~ se m m 
% 4» 48 4€ V 6S 68S £N 48 съ 49S 6% 0% дъ съ съ дъ дъ Фъ 78 048 0€ дъ Фъ Фъ Фъ #9 дъ съ дъ EE EL FH ED sR eH ER ER Ed 48 дъ дъ 48S 48 


7, ра дара ду даа ра mm p» МАЈА ара а 
дъ Фъ дъ Фу #у 27 4> >: Фу 


ET = 0.1 
праз 

ЕТ = Е 
EIE 


pog e Ce oDFOO/C(MM А) 
(002 -0B*OO0*RJ/A) - CCI = (Е **2))/ET)SDCOSC(TA))-* FR - 
ГЕТЕ си и *AW*ET)*(1- R/CÀ*(C1-(E**2))))"*(SINCTA)*FS)- 
(T * DMM) 
FETURN 
END 
E a es Се АА 


еее еее еее ее 962727672267 


* + CALCULATE THE и ОКВТТАТ, ELEMENTS 
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SUBROUTINE NSMA(A,DT,DA,NEWA) 
THIS SUBROUTINE CALCULATES THE NEW SEMI-MAJOR AXIS 


DOUBLE PRECISION DA,DT,A,NEWA 


NEWA = A + DA*DT 


RETURN 
END 


еее: 


SURSOUDENE NECOCE,DT,DE, NEWE) 
es SUBROUTINE CALCULATES THE NEW ECCENTRICITY 


РЕКЕ PRECISION DE,DT,E,NEWE 
ME E -DE*DT 

RETURN 

ENT) 


TVET OTE TS Te Te eRe TPO EE TES 


SUBROUTINE NINCL(1I,DT,DI,NEWI) 

Mts SUBROUTINE CALCULATES THE NEW INCLINATION 

ШИРЕЕР РКЕСТЪТОЖ ВТ ОТ, Г, към т 

ЕТ = Те DI*DT 

ТЕСЕК 

ЕХО 

ое ae oR CIC Cae te 


ас У У УУ У У У У У У У 


SUBROUTINE NASNOD( LAN , DT ,DLAN , NEWLAN) 


ORB18420 
ORB18430 
ORB18440 
ORB18450 
ORB18460 
ORB18470 
ORB18480 
ORB18490 
ORB18500 
ORB18510 
ORB18520 
ORB18530 
ORB18540 
ORB18550 
ORB18560 
ORB18570 
ORB18580 
ORB18590 
ORB18600 
ORB18610 
ORB18620 
ORB18630 
ORB18640 
ORB18650 
ORB18660 
ORB18670 
ORB18680 
ORB18690 
ORB18700 
ORB18710 
ORB18720 
ORB18730 
ORB18740 
ОКВ18750 
ОКВ18760 
ОКВ18770 
ОКВ18780 
ORB18790 
ORB18800 
ORB18810 
ORB18820 
ORB18830 
ORB18840 
ORB18850 
ORB18860 
ORB18870 
ORB18880 
ORB18890 
ORB18900 


THIS SUBROUTINE CALCULATES THE NEW LONGITUDE OF THE ASCENDING NODEORB18910 


DOCSDESERBOISION РТАМ,рТ, ТАМ, МЕКА 


NEWLAN = LAN + DLAN*DT 
Ecos 
ES 


с 
NS) 


ORB18920 
ORB18930 
ORB18940 
ORB18950 
ORE18960 
ORB 13970 


453696329227 


3% 
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SUBROUTINE NARPERCAP,DT,DAP,NEWAP) 
THIS SUBROUTINE CALCULATES THE NEW ARGUMENT OF PERIGEE 


DOUBLE PRECISION DAP,DT,AP,NEWAP 
NEWAP = AP + DAP*DT 

RETURN 

END 


veoesezeyesevyes ee vede eyevevevevedevede dev edeses ee yeyeyedyeve vere Fede re yes ev eye eee ye ye ye ye ve sede ve 


SUBROUTINE NMNAN(MA,MM,DT,TF,DMA,PI) 
THIS SUBROUTINE CALCULATES THE NEW MEAN Anomaly 


DOUBLE PRECISION DMM ВК В5, ВМА ОТО ДВЕТЕ КИ ee ГЕН 


JA = УМУСТЕ) F DIAD 


TEMA -CGT А УБ TOS EN 
МА = MA- (2"PT) 

ЕЗІЛЕ 

RETORN 

END 


УУ ИУ о reeek ee aeree 
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SUBROUTINE NMNMOCMM, DMM, DT, NEWMM) 

THIS SUBROUTINE CALCULATE ТНЕ МЕИ ИЕА ИОД: 
DOUBLE PRECISION ОНА О РА РЕМСА 
КЕЧИН = М + римкрт 

RETURN 


EN 
АЕ 


-- ec'anlen afa n'a n РТА РРА РА РРА РРА РРА РТА РТА РТА РТА РУМ РТА РТА РУТЕ РРА РУК РУНО nls ula u'a afa nlanla alanla aea, Sad, моја јаје ја 5% edu ula uj, sese o uiua udo uon РАЊАВАЊА 
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SUBROUTINE NEA(MA,E,EA) 

THIS SUBROUTINE CALCULATES THE NEW ECCENTRIC ANOMOLY BY USING 
NEWTONS METHOD OF ROOT FINDING 

DOUBLE PRECISION EAN,MAN,MA,E,EA,DIFF 

LET (EA) EQUAL (MA) FOR INITIAL GUESS AT ROOT 

EA = MA 

EAN = EA + (MA - EA +E*DSIN(EA))/(1.0D+00 - E*DCOS(EA)) 
MAN = EAN - E*SIN(EAN) 


CHECRSVIFFERESCE MCDTEF) 
DIFF = АВЗСМА -МАМ) 
ЕА = ЕАМ 


CONTISCE ТО ТУТЕКАТЕ UNTIL DIFFERENGE ITS МЕБ ВОТЬ Е 
ТВ ОГ, 0. 90009000001) ЕН 
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ORB18980 
ORB18990 
ORB19000 
ОВО 
ORB19020 
ORB19030 
ORB19040 
ORB19050 
ORB19060 
ORB19070 
ORB19080 
ORB19090 
ORB19100 
ОКВ19110 
CRP 13 E20 
ORB19130 
ORB19140 
ОВВ1 915 
ОКВ19160 
ORB19170 
ORB19180 
ORB19190 
ORB19200 
ORB T2210 
ORB 19220 
ORB19230 
ORB19240 
ОКВ19250 
ОКВ19260 
ОКВ19270 
ОКВ19280 
ORB19290 
ORB19300 
ORB 193 KO 
ORB 19320 
ОКВ 19330 
ОКВ 19340 
ORB19350 
ORB19360 
ORB19370 
ORB19380 
ORB19390 
ORB19400 
ORB19410 
ORB19420 
ORB19430 
ORB19440 
ORB19450 
ORB19460 
ORB19470 
ORB19480 
ORB19490 
ORB19500 
ORB19510 
ORB19520 
ОКВ 


в А C EASESECSDSIDNCEA)9Z( 1.00500 s. E*DCOSCEA)) ORB19540 

MAN. 2 EAN"- E*DSIN(EAN) ШЕР 9 50 

EA = EAN ORB19560 

DIFF = ABS(MA = MAN} 0819570 

GOTO 200 ORDISSSO 

ВЕ ОКВ 19590 

EA = EAN ORB19600 
RETURN ORB19610 

ESD ORB19620 
ORB19630 

Aireen e erer ere ere e еее ieee kek kekeke rA eae eeke e e eeek ORB19640 
ORB19650 

SUBROUTINE NTACEA,E,TA,PI) ORB19660 

72 THIS SUBROUTINE CALCULATES THE NEW TRUE Anomaly ORB19670 
ORB19680 

DOUBLE PRECISION EA,E,TA,PI ORB19690 
ORB19700 

Та - DACOSCCE - DCOS(EA))/(E*DCOS(EA) - 1.0D*00)) ORS 197 10 
PEREA. GI- PI) THEN ОКВ19720 

= (2*PI) - TA ORB19730 

ENDIE ORB19740 
RETURN ORB19750 

END ORB19760 
ORB19770 

Поа ти и СС СС Со о У Со genes esesesevesevese sese oe ve dyes esyeveseoeveve DCIS IC DS IS IIS IT TE У Је ORB19780 
ORB19790 

SUBROUTINE NANGHMO(H,DT,DH,NEWH) ORB19800 

* iis SUBROUTINE CALCULATES THE NEW ANGULAR MOMENTUM ORB19810 
ORB19820 

ESUBBESOPRECESION DH,DT,H;NEWH OKBI9S830 
ORB19840 

р = НН = БЕРТ ORB19850 
RETURN ORB19860 

END ORB19870 
ORB19880 

ке С СГ С ССИ ГГ С ТС СС С С С СИ С СК СК С С У СУХИ ЋИР ИУУ ORB19890 
ORB19900 

SeeROUlINE INITSUM( TFEA,TESU,TEMO,TPDRA,TDI,IDA,TDE,TDMM,TDMA, ORB 19910 

+ TDLAN, TDH TDAP) ORB19920 

* THIS SUBROUTINE INITIALIZES THe sUMS OF £ORCES AND ELEMENT CHANGESORB 19930 
ORB19940 

ПАО 6 ОКВ19950 

TFSU = 0.0 ORB19960 

TFMO = 0.0 ORB19970 

TrDRA = 0.0 ORB19980 

ТО + 0.0 ОКВ 19990 

TDA = 0.0 ORB20000 

ПОЕТ = 070 ОКВ20010 

тр 0-0 ORB20020 

TDMA = 0.0 ORB20030 

TDLAN = 0.0 ORB20040 

IE 0 ORB20050 

TDAP = 0.0 ORB20060 

eee UN ORB20070 

Ваш ORB20080 
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sese yes eese vevede desee veg esee ev edesev ee deve desedese eee еее еее еее еее: 


T CALCULATE me _ NEN POSITION ea ыы, а 


БАНИ ООН сУ 3ezeveses esee esee ese ves eses esee deveveeve esee eese esee eo 


SUBROUTINE NPOS(RI,RJ,RK,R,LAN,AP,INC, TA,A,E) 


THIS SUBROUTINE CALCULATES THE NEW POSITION VECTOR 

DOUBLE PRECISION XW,YW,ZW,INC,RI,RJ,RK,R,LAN,AP,TA,A,E 
“e CALCULATE POSITION VECTOR IN 'PQW' SYSTEM 

= (A*(1 - (E**2)))/(1 + E*DCOS(TA)) 

XW = R*DCOS(TA) 

W = R*DSIN(TA) 

ZW = 0 
* TRANSFORM POSITION TO 'IJK' SYSTEM 


CALL PQWIJK( LAN, AP, INC, AW, eW, ZW RR 
a DSQRT((RI**2) + (ЕЛ2) - (КК 2)) 
RETORN 

END 


А eye yeveye sese sees e eyes kekke rke kkk kekere kekk 


SUBROUTINE NVEL(E, P, TA ВА, АР, ME heavenly) 
THIS SUBROUTINE CALCULATES THE NEW VELOCITY VECTOR 


DOUBLE PRECISION- INC VE VO NW MO Er р т де КО 


CALCULATE, VELOCITY IN § POW SYSTEM 
VP DSQRT( MU/P)*( -DSIN(TA)) 
DSORTCHU/POSUCESE ODOOSCTAD) 

D DEO 


< 
Su 
пи | 


* TRANSFORM VELOCITY INTO IJK’ SYSTEM 
CALL POWIJK(LAN AP INCS VP в) 
= DSORT( Ут 27) + ( VJ****2) қ VK**2 ) 
RETURN 
END 


Vevesedevedevese ses ezevyesezed ese eyeevevyedezese eyes e dedede ye Jeedyeve ede yes e deve sedes УУ У УУ У У У По ev ede deve deve 


VELOCITY CHANGE 
ТЕ КЕСЕК КЕСЕК СЕЕЕСЕККЕСТЕСЕЕТЕСІІ 25907692522 6925565696927290922656909996909699с972929256963692596969055909696969696969296969696929 292 


SUBROUTINE CHGVEL(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R, 
+ VI,VJ,VK,V,MU,PI,H,A,E,N,TA,P,MM,MA,EA, 
+  TF,T,NUM,RIRAY,RJRAY, RKRAY, RARAY, TARAY, AINRAY, APRAY, TIMRAY, 
+  TT,EI,EJ,EX,LP,HI,HJ,IOPT1,TFEA,TFSU,TFMO,TFDRA, 
+  TDI,TDA,TDE,TDMM,TDMA, TDLAN, TDH, TDAP) 
> THIS SUBROUTINE CALCULATE VELOCITY CHANGES 


* THE FOLLOWING SUBROUTINES ARE CALLED: 
TACHG = RETURNS TRUE ANOMALY FOR VELOCITY CHANGE LOCATION (CHTA) 


AND AN INDICATOR OF LOCATION (ITA) 
CALCEL = CALCULATE Orbital ELEMENTS 
CAPRET = CALCULATE GNEERTURBED ORBIT 


ORB20090 
ORB20100 
ORB20110 
ORB20120 


ORB20130 | 


ORB20140 
ORB20150 


ORB20160 © 


ORB20170 
ORB20180 
ORB20190 
ORB20200 
ORB20210 
ORB20220 
ORB20230 
ORB20240 
ORB20250 
ORB20260 
ORB20270 
ORB20280 
ORB20290 
ORB20300 
ORB20310 
ORB20320 
ORB20330 
ORB20340 
ORB20350 
ORB20360 
ORB20370 
ORB20380 
ORB20390 
ORB20400 
ORB20410 
ORB20420 
ORB20430 
ORB20440 
ORB20450 
ORB20460 
ORB20470 
ORB20480 
ORB20490 
ORB20500 
ORB20510 
ORB20520 
ORB20530 
ORB20540 
ORB20550 
ORB20560 
ORB20570 
ORB20580 
ORB20590 
ORB20600 
ORB20610 
ORB20620 
ORB20630 
ORB20640 


г. 
Сл 


а 
+ 
ще 
+ 
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NPOS = CALCULATE NEW POSITION 

NVEL = CALCULATE NEW VELOCITY 

STORE = STORE POSITION AND ELEMENTS IN ARRAYS 

ENERGY = EXERGY OF SATELLITE 

ECC = ECCENTRICITY 

SHAXIS — SEMI-MAJOR AXIS 

COUBLE PRECISTON T. DI PER AL EAN AP ERIS RI RKR; VI VJ, VR; V, 

MORE RA EN TASR M MAEA TIE; IT; 

ESPENMISTNENVIONEWVKSNENVSVMAR,CHTASEISEJ,;EK,LP,HIJHJ,VCIR, 
DI,DE,DA,DMM,DMA,DLAN,DH,DAP,NEWEI,NEWEJ,NEWEK,NEWE , NEWENR, 
NEWA,NEWRP, RE 


DIMENSION RARAY(500),TARAY(500) ,RIRAY(500) ,RJRAY( 500), 
RKRAY( 500) , AINRAY( 500) , APRAY(500) , TIMRAY( 500) 


CHARACTER*1,YORN, PYORN 
КЕ = 6. 37820+03 


PROMPT THE USER FOR THE VELOCITY Change LOCATION 
CALL TACNG(PI,CHTA, ITA) 


ИРЕ СООЦГОТЕК TO ОКЕ ТТМЕ STEP 
T = DT 


ROTATE TO THE VELOCITY CHANGE LOCATION 
THIS IS IDENTICAL TO THE Unperturbed ORBIT WITH THE EXCEPTION 
THAT A COMPLETE ORBIT IS NOT CALCULATED 
PRINT*, ROTATE TO VELOCITY CHANGE LOCATION’ 
ШІ С(ТТА ЕП 2) ОБ. (ІТА. ЕО. 3)) ТНЕМ 
PRINT*,'BEFORE TA =',ТА 
IF (TA .GT. 6.21) THEN 
A= TAS OPI) 
ENDIF 
ial Gober) ASADA TA. LIT. CHTAJ) THEN 
PRINT*,'TA =',TA 
NUM = NUM + 1 
TT = TT + DT 
CALL NEWELT(MM,MA,E,EA,TA,TF,DT,PI, PER) 
CALI NPOS(RI RI, RK, R LAN AP, 1, TALASE) 
CALL AVELL E P TAS LAN APIT VI VJ VE, V MU) 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY, RARAY, 


t TARAY,NUM,I,AP,AINRAY,APRAY,TT, TIMRAY) 


T = T + DT 
Сото 230 
ENDIF 
ПР Do GEM PER) IMEN 
DEDE PER 
ЕМЕ 
ENDIF 


PRINT ESCAPE VELOCITY AND CIRCULAR VELOCITY FOR Reference 
в а Бро СЕ) 

РЕР А. = TA 

PRINT*,'THIS SHOULD BE THE DESIRED RADIUS RP OR RA’ 


ORB20650 
ORB20660 
ORB20670 
ORB20680 
ORB20690 
ORB20700 
ORB20710 
ORB20720 
ORB20730 
ORB20740 
ORB20750 
ORB20760 
ORB20770 
ORB20780 
ORB20790 
ORB20800 
ORB20810 
ORB20820 
ORB20830 
ORB20840 
ORB20850 
ORB20860 
ORB20870 
ORB20880 
ORB20890 
ORB20900 
ORB20910 
ORB20920 
ORB20930 
ORB20940 
ORB20950 
ORB20960 
ОЕВ20970 
ORB20980 
ORB20990 
ORB21000 
ORB21010 
ORB21020 
ORBZ 1030 
ORB21040 
ORB21050 
ORB21060 
ORB21070 
ORB21080 
ORB21090 
ORB21100 
ORB2TTIO 
ORB21120 
ОВО 
ORB21140 
ORB2 T130 
ОБВ2 1160 
ORB21170 
ORB21180 
ОВА ИО 
ОК 2200 


260 PRINT, 'RADIUS 


=! NS 
PRINT... VELOC] iia E 
VMAX = DSORT( 2. 0 (1U EERI 
PRINT*,'MAX VELOCITY AT THIS RADIUS То 
VCIR = DSQRT(MU/R) 


PRINT*, CIRCULAR VELOCITY AT THIS ҚАР Б OUR 


PROMPT USER TO CHANGE VELOCITY IN ORBITAL PLANE 


PRINT*,'DO YOU WANT TO CHANGE THE VELOCITY IN THE ORBITAL PLANE?' 


PRINT*,'ENTER "Y" OR "N" :' 
READ* , PYORN 

PRINT* , PYORN 

ТЕ CPYORN ЕО ЭТНЕ 


PRINT*, 'GIVE THE TOTAL CHANGE IN VELOCITY, I.E. 5.0 KM.' 
PRINT*,'THE PROGRAM WILL FIGURE OUT THE FINAL VELOCITY VECTOR’ 


PRINT*,' ENTER VELOCITY CHANGE: ' 
READ* , CHGV 
PRINT* , CHGV 


CALCULATE NEW VELOCITY FOR CHANGE IN THE ORBITAL PLANE 


NEWVI = VI + (CHEV МЕТИ) 
NEWVJ = VJ + (CHGV * VJ / V) 
NEWVK = VK + (CHGV * УК / У) 


Velocity CHANGE OUT OF ORBITAL PLANE 
ELSEIF (PYORN ЕО 'N') THEN 
PRINT*,' ENTER THE NEW VELOCITY VECTOR: ' 
PRINT*,' ENTER THE NEW VI' 
READ* , NEWVI 
PRINT* ,NEWVI 
PRINT*,' ENTER THE NEW VJ' 
READ* , NEWVJ 
PRINT*,NEWVJ 
PRINT*,' ENTER THE NEW VK' 
READ*,NEWVK 
PRINT*,NEWVK 


NUM S] 
ГТ VOS 

ELSE 
CALL EXCMS('CLRSCRN') 
GOTO 260 

ENDIF 


PRINT NEW VELOCITY FOR USER TO CHECK 
NEWV - DSQRT((NEWVI?^*2) + (КЕКУЛ 2) + (NEWVK***2)) 
РЕТМТ“, NEW VI =' ,NEWVI 
PRINT*,'NEW VJ z' ,NEWVJ 
PRINT*,'NEW VK = ,NEWVK 
PRINT, 'NEW Vo =' ,NEWV 
PRINT*,' ARE THESE VALUES THE ONES YOU WANT?' 
PRIMIS ЗЕР "y ОКО Neen 
READ* , YORN 
PRINT* , YORN 
ТТК (ҮСЕМ SEC "UN aes 
СӘТІ ЕШСМ5( "СІЗБЕН 
GOTO 260 


64 


ОКВ 1216 
ORB21220 
ORB21230 
ORB21240 
ORB21250 
ORB21260 
ORB21270 
ORB21280 
ORB21290 
ORB21300 
ORBZISTD 


'ORB2 1320 


ORB2 1550 
ORB21340 
ORB21350 
ORB21360 
ORB 21370 
ORB2 1380 
ORB21390 
ORB21400 
ORB21410 
ORB21420 
ORB21430 
ORB21440 
ORB21450 
ORB21460 
ORB21470 
ORB21480 
ORB21490 
ORB21500 
ORB21510 
ORB21520 
ORB21530 
ORB21540 
ORB21550 
ORB21560 
ORB21570 
ORB21580 
ORB215390 
ORB21600 
ORB21610 
ORB21620 
ORB21630 
ORB21640 
ORB21650 
ORB21660 
ORB21670 
ORB21680 
ORB21690 
ORB21700 
ОКВ2 710 
ОКВ21720 
ОКВ (20 
ОКВ2 1740 
ОКБ ЗИ 
ОКВ2 1760 


DA 


4 


ate 
» 


+ + + 


РАТЕ 


ECR FOR VALID VELOCITY 

ТЕРЕНИ СІ. VHAN) THEN 
PRINT*,' YOUR VELOCITY IS TO GREAT !!' 
GOTO 260 

ENDIF 


Calculate PERIGEE RADIUS TO SEE IF SATELLITE WILL IMPACT EARTH 
CALL ENERGY(NEWV,R,MU,NEWENR) 
CALL ECC(RI,RJ,RK,R,NEWVI,NEWVJ,NEWVK,NEWV ,NEWEI,NEWEJ , NEWEK, 
NEVE , MU) 
CALL SMAXIS(MU,NEWENR,NEWA) 
NEWRP 2 NEWA*(1.0 - NEWE) 
IF (NEWRP . LE. RE) THEN 
PRINT*,'YOUR VELOCITY AT THIS POINT IS TO SMALL!!!' 
PRINT*,'THE SATELLITE WILL IMPACT THE EARTH!!' 
PRINT*,'THE SATELLITES RADIUS OF PERIGEE WOULD BE ',NEWRP 
PRINT*,'A NEW VELOCITY WILL HAVE TO BE ENTERED!! ' 
GOTO 260 
ENDIF 


ACCEPT NEW VELOCITY 
VI NEWVI 

УЈ = МЕКУЈ 

VK = NEWVK 

\ = NEWV 


CALCULATE NEW ELEMENT WITH NEW VELOCITY AND SET TIME STEP 

BESUNCALGELC RI,RJ,RK,R,VI,VJ,VK,V,EI,EJ,EK,E, A, I, LAN, LP , TA, 
PER,EA,MA,AP,AL, TF, P, PI, MU, MM, N, H, HI,HZ) 

РТ = РЕК/50.0 

T = DT 


THE FOUR Different CASES OF VELOCITY CHANGES FOLLOWS: 


VELOCITY CHANGE AT PERIGEE, AND NEWV > V Circular 
IF(( ITA. EQ. 1). AND. (NEWV. GT. VCIR) ) THEN 
CALL UNPRET(DT,PER,AL,LAN,AP, I,RI,RJ,RK, R, VI, VJ, VK, V, 
MU,PI,H,A,E,N, TA, P, MM, 
MA,EA,TF,T,NUM, RIRAY , RJRAY , RKRAY , RARAY, 
TARAY , AINRAY , APRAY, TIMRAY,TT) 


Change VELOCITY AT PERIGEE, AND NEWV «- V CIRCULAR 
APOGEE AND PERIGEE SWAP 
ШЕБЕР COllA. EQ. 1). AND. (NEWY. LE. VCIR) )THEN 


CLEAR PREVIOUS PLOTS 
NI 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY , RARAY , TARAY, 


T NUM,I,AP,AINRAY,APRAY,TT, TIMRAY) 


FER? 


$ 


Sete ELLIE TO NEW PERIGEE- -SONLY 4AJHALERORBIT 
ОЕЕО Е БР, ВЕ УЕ, КУ, 


+ СЦЕ ESSA ER, P T, 


ORBZ17 70 
ORB21780 
OR B27 90 
ORB21800 
ORB21810 
ORB21820 
ORB21830 
ORB21840 
ORB21850 
ORB21860 
ORB21870 
ORB21880 
ORB21890 
ORB21900 
ORB21910 
ORB21920 
ORB21930 
ORB21940 
ORB21950 
ORB21960 
ORB21970 
OREZ 1280 
ORB21990 
ORB22000 
ORB22010 
ORB22020 
ORB22030 
ORB22040 
ORB22050 
ORB22060 
ORB22070 
ORB22080 
ORB22090 
ORB22100 
ОР ТО 
ORB22120 
ORB22130 
ORB22140 
ORB22150 
ORB22160 
ORB22170 
ORB22180 
ORB22190 
ORB22200 
ORB22210 
ORB22220 
ORB22230 
ORB22240 
OREZ2250 
ORB22260 
ORB22270 
ORB22280 
ORB22290 
ORB 22300 
ОКВ2 2310 
ОЕ О 


+ + + 


MA,EA,TF,1T,NUM,RIRAY , RJRAY , RKRAY, RARAY, 
TARAY , AINRAY , APRAY , TIMRAY , TT) 


RESET TIME COUNXTERODESONESBRUESSTER 
T = DT 


CALCULATE COMPLETE NEXT ORBIT 
CALL UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ,VK.V, 
MU,PI,H,A,E,N, TA, P,MM, 
MA,EA,TF,T,NUM,RIRAY, RJRAY , RKRAY , RARAY, 
TARAY , AINRAY , APRAY , TIMRAY , TT) 


CHANGE VELOCITY AT APOGEE, AND NEW V < V CIRCULAR 
ELSEIF (CITA. EQ. 2) . AND. CNEWV „БТ. УЕ! 


+ + + 


TX "PER? Z 


ВИНЕСЕ В 
CALL УМРАЕТЕОТ, РЕВ, АБ, БАМ, АР, ВТ ВТ ВЕ И 
МОЭРТН,А, Б № 
MA,EA, TF,T,NUM,RIRAY,RJRAY,RKRAY,RARAY, 
TARAY , AINRAY , APRAY , TIMRAY , TT) 


CHANGE VELOCITY AT Apogee, AND NEWV >= V CIRCULAR 
OR AT ANY OTHER TRUE Anomaly 
ELSEIF ( (CITA. EQ. 2). ANDZOCNEWV. бЕГУСІК)) ОКЕ ОТТА. СОИ 


+ + 


+ + + 


ТЕ (ТА бо КНП 
ТА = TA - (2*PI) 
ENDIF 


CLEAR PREVIOUS ORBITS AND STEP SATELLITE TO ЗБ ВЕ 
1- IF 
ACT =l 
CALL STORE(RI,RJ,RK,R,TA,RIRAY,RJRAY,RKRAY,RARAY,TARAY, 


NUM, I, AP, AINRAY, APRAY , TT, TIMRAY) 
CALL UNPRET( DT, PER, AL, LAN, AP, I,RI,RJ,RK,R, VI, VJ,VK,V, 
MU,PI,H,A,E,N, TA, P,MM, 
MA,EA,TF,T,NUM, RIRAY, RJRAY , RKRAY , RARAY, 
TARAY , AINRAY , APRAY , TIMRAY , TT) 
IF (TF .GE. PER) THEN 
TF = TF - PER 
ENDIF 
CALCULATE COMPLETE NEXT ORBIT 
Т = DT 
CALL UNPRET(DT,PER,AL,LAN,AP,I,RI,RJ,RK,R,VI,VJ, VK, V, 
MU,PI,H,A,E,N, TA, P,MM, 
MA,EA,TF,T,NUM,RIRAY,RJRAY , RKRAY, RARAY, 
TARAY , AINRAY , APRAY, TIMRAY,TT) 


8-219 9. 9 —* е ео 
--.-.-................”-.. 


ORB22330 
ORB22340 
ORB22350 
ORB22360 
ORB22370 
ORB22380 
ORE2 2390 
ORB22400 
ORB22410 
ORB22420 
ORB22430 
ORB22440 
ORB22450 
ORB22460 
ORB22470 
ORB22480 
ORB22490 
ORB22500 
ORB225 80 
ORB22520 
ORB22530 
ORB22540 
ORB22550 
ORB22560 
ORB22570 
ORB22580 
ORB22590 
ORB22600 
ORB22610 
ORB22620 
ORB22630 
ORB22640 
ORB22650 
ORB22660 
ORB22670 
ORB22680 
ORB22690 
ORB22700 
ORB22710 
ORB22720 
ORB22730 
ORB22740 
ORB22750 
ORB22760 
ORB22770 
ORB22780 
ORB22790 
ORB22800 
ORB22810 
ORB22820 
ORB22830 
ORB22840 
ORB22850 
ORB22860 
ORB22870 
ORB22880 


uo n'a eon n asosoy 
4% 4% 428 98 4S 4S 6 S 
ыы OUTPUT PLOTS 
ale m - m ela o's m m a's a's ale 
wee ee ar ee ed 4% 48 оъ Ф 


SUBROUTINE TACNG(PI,CHTA,ITA) 
THIS SUBROUTINE Asks THE USER FOR VELOCITY CHANGE LOCATION 


POCE CE PRECISION CHTA,PI 


CALL EXCNSC CLRSCRN") 

PRINT*,'WHERE DO YOU WANT TO CHANGE THE VELOCITY? ' 
PRINT*,' 1. AT CURRENT PERIGEE' 

PINT, 2. AT CURRENT Apogee’ 

PRINT*,' 3. AT À SPECIFIC TRUE Anomaly' 

hide, ENTER 1, 2 ОВ 3” 

READ* , ITA 

PRINT*,ITA 


SET TRUE ANOMALY CHANGE LOCATION (CHTA) TO DESIRED LOCATION 
Mella, EQ. 1) THEN 


CHTA = 0.0 

EDIF 

ШЕСТА СВО. 2) THEN 
СНТА = PI 

ENDIF 


ПИ ШТА ЕО 3) THEN 
PRINT*,'AT WHAT TRUE ANOMALY DO YOU WANT TO CHANGE THE' 
ГЕ VELOCITY? ' 
PRINT*,' ENTER TRUE ANOMALY IN DEGREES' 
READ* , CHTA 
PRINT*,CHTA 
ЛЕС ons PI 180 
ENDIF 
RETURN 
END 


emu au au Dn. 


Мама УУ 


а Vesesegeseseveveseveveses ese ves esee ves eyeveves ees ese eye des eee eye селе ее ттт 


=== 
27... 


SUBROUTINE PLOTS(RIRAY,RJRAY,RKRAY,RARAY,TARAY ,NUM,PI,INC,LP,A, 
E ENDE AINKRAYJAPRAY /IDIUBRAY DEEASTESUSTEMOSTEDRA; 
+ РОБА Oly Dente Ge uA, IDLA IDH; IDAR; 
4 


META LASN HTAR R, Y 


BEEN SSSUBEOUPINE ASKS THE USER FOR THE TYPE OF OUTPUT THAT 15 
ШЕ ІКЕЙ EERIFOCAL, GROUND TRACK OR TO SKIP THE PLOT. 


THE FOLLOWING SUBROUTINES ARE CALLED: 


PERIF = PLOT PERIFOCAL ORBIT 

GRTRK = PLOT GROUND TRACK 

РАТЕ = DISPLAYS DATA ON PLOT 

Ше ос ubusspla TO TEC 618 OUIPUT 
ENDPL = END THIS DISSPLA PLOT 


REFER TO DISSPLA USER'S MANUAL FOR EXPLANATION OF DISSPLA 
SUBROUTINES 


Dec MESENECISIONOPISAE INCDBE,TE,FER;MM HA, LAN, H, AP,R,V 


ORB22890 
ORB22900 
ORB22910 
ORB22920 
ORE2Z 2950 
ORB22940 
ORBZ2950 
ORB22960 
ORB22970 
ORB22980 
ORB22990 
ORB23000 
ORB23010 
ORB23020 
ORB23030 
ORB23040 
ORB23050 
ORB23060 
ORB23070 
ORB23080 
ORB23090 
ORB23100 
ORB23110 
ORB23120 
ORB23130 
ORB23140 
ORB23150 
ORB23160 
ORB23170 
ОКВ 3180 
ORB23190 
ORB23200 
ORB23210 
ORB2 3220 
ORB23230 
ORB23240 
Ве Ао 
ORB23260 
ORB23270 
ORB23280 
ORB 23290 
ORB23300 
ORB23310 
ORB23320 
ORB23330 
ORB23340 
ORB23350 
ORB23360 
ORB23370 
ORB23380 
ОКБ? 590 
ОКВ2 3400 
ORB23410 
ORB23420 
ORB23430 


DIMENSION RIRAY(500) ,RJRAY( 500) , RKRAY(500) ,RARAY( 500) , TARAY( 500), ORB23440 


AINRAY( 500) ,APRAY( 500) , TIMRAY( 500) 
CHARACTER*1, YORN 


CATEEESCHUISC CLRSORSEM 


* CALCULATE SINGLE PRECISION VARIABLES 


(2 


340 


SPI = SNGLC(PI) 


SA = SNGL(A) 

SE - SNGL(E) 
SINC = SNGLCINC) 
SLP = SNGLCLP) 
STF = SNGLCTF) 
SPER = SNGL( PER) 
SMM = SNGL(MM) 
SMA = SNGL(MA) 


SLAN = SNGLC LAN) 
SH = SNGL(CH) 

SAP = SNGLCAP) 
SV = SNGLC(V) 

SR = SNGL(R) 


PROMPT USER FOR DISECA? PITRE 
PRINT*,' WHAT TY Ree Display 15 DESIRED 


PRINT*, 1. PERIFOCAL 
PRINT*,' 2. GROUND TRACK' 
PRINT*,' 3. SKIP PLON 
PRINT ENTER а ан 
READ” , INPUT 

РЕТ ОКО РИ 

ЕОКМАТ( 14) 


САБЫ ТЕР 18 


CALL APPROPRIATE PLOT 
IF (INPUT EO. 1) THER 

CALL PERIF(RARAY,TARAY,NUM,SPI,SINC,SLP,SA,SE) 
ELSEIF (INPUT .EQ. 2) THEN 

CALL GRTRK(AINRAY,APRAY,TARAY,STF ,NUM, TIMRAY) 
ELSEIF (INPUT .EQ. 3) THEN 

GOTO 360 
ELSE 

РЕГ, INVALID En 

GOTO 340 
ENDIF 


DISPLAY DATA 
САЪЬ БАТАСЗТМС,5А,5Е,ТЕЕА, ТЕЗУ, ТЕМО, ТЕЮБА, РЕВ. БРУТ, ОДЕТ 


+ ТЬМУ, ТЬМА , ТЬБАМ, ТВН, ТРАВА, ОЛИО SIS T ПОН 


CALLTENDPLC(O) 


PROMPT USER IF ANOTHER DISPLAY TYPE IS DESIRED 
PRINT*,'WOULD YOU LIKE ANOTHER PLOT USING THE SAME ORBITAL' 
PRINT*,'PARAMETERS AND DATA: ' 

RINT ENTER т Ge Па E. 
READ* , YORN 

PRINT*,YORN 


ORB23450 
ORB23460 
ORB23470 
ORB23480 
ORB23490 
ORB23500 
ORB23510 
ORB23520 
ORB23530 
ORB23540 
08823558 
ORB23560 
ORB23570 
ORB23580 
ORB23590 
ORB23600 
ORB23610 
ORB23620 
ORB23630 
ORB23640 
ORB23650 
ORB23660 
ORB23670 
ORB23680 
ORB23690 
ORB23700 
ORB237 1¢ 
ORB23720 
ORB23730 
ORB23740 
ORB23750 
ORB23760 
ORB23770 
ORB23780 
ORBZ37 30 
ORB23800 
ORB23810 
ORB23820 
ORB23830 
ORB23840 
ORB23850 
ORB23860 
ORB23870 
ORB23880 
ORB23890 
ORB23900 
ORB23910 
ORB2 3920 
ORB23930 
ORB23940 
ОКВ23950 
ORB23960 
ORB23970 
ORB23980 
окваз а о 
ORB24000 


| ЖУО ЕО У’) ТЕМ 
GOTO 340 
ENDIF 
360 RETURN 
END 


424 280562 


EDPEROUITSE PERIFP(RARAY, TARAY,NUM,PI,INC,LP,A,E) 
* О СОЕКОО МЕ ОРОТ ООТ НЕ KESULIS Of THE PROGRAM USING THE 
=~ Wise LAY FEATURE ON THE MAIN FRAME. 
E REFER TO DISSPLA USERS GUIDE FOR EXPLANATION OF DISSPLA 
* OCBROUTINES. 


PEAL INC, LP 
DIMENSION TARAY(500),RARAY(500),RIRAY(500),RJRAY(500),RKRAY(500) 


I = 1 


s SET SCALE OF AXIS 
RSTEP = (A*(1+E)) / 3 
CALL TEK618 
CALL RESET(3HALL) 
SCMPLX 
БАРІ PHYSOR(1.25,4.) 
CALL AREA2D(6. ,6. ) 
CALL MESSAG('PERIFOCAL COORDINATE SYSTEMS’ ,100,1.0,6. 5) 
в В АЕ W ,2) 
РАТЕ YNAME( YW’ ,2 
CALL XAXANG( 90. 0) 
GALL YAXANG(O. 0) 
INTAXS 
POLAR ЈА RSTEP.3. ,3. ) 
POLY3 
CALL NOCHEK 
CALL CURVE(TARAY,RARAY ,NUM, 1) 
BALL СОМРЕХ 
CALL HEIGHT(. 2) 
CALL RESET( 'COMPLEX' ) 
CALL RESET( HEIGHT!) 
CALL ENDGR(O) 


* Display EARTH PLOT 
CALL EARTH1(A,E, INC,PI,RSTEP) 
RETURN 
END 


4.................../.......... Jc m m о с PEE PES PPS PPS PER PPS PM 
PL ERR ER OD CN es ER ES Ed GL GU EDL EL EN EL CS PR EX ER EN CR ED CL 048 ду дъ Фу Фу гу 


е а а ја а о ова е ъв но пао оъ шј йо йо го је во ъв =з 
6% 9% ex ex ex ev eb vC vs EN 4% 4% дъ Фа Фу дъ съ 48 Фу Фа Фъ гу 


EEPRODDBENESEARTHTCA,E,INCSPIRSTEP) 
х THIS SUBROUTINE PLOTS A VIEW OF THE WORLD, LOOKING DOWN THE 'Z' 


з Дао, РЪЛСЕВ ОТНЕ ОКТОТУ — НЕ Latitude IS FIXED, BUT THE 
3 LONGITUDE VARIES WITH THE INCLINATION. 
e REFER TO DISSPLA USER'S MANUAL FOR EXPLANATION OF DISSPLA 


ЗОВОО ЕЕ 
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ORB24010 
ORB24020 
ORB24030 
ORB24040 
ORB24050 
ORB24060 
ORB24070 
ORB24080 
ORB24090 
ORB24100 
ORB24110 
ORB24120 
ORB24130 
ORB24140 
ORB24150 
ORB24160 
ORB24170 
ORB24180 
ORB24190 
ORB24200 
ORB24210 
ORB24220 
ORB24230 
ORB24240 
ORB24250 
ORB24260 
ORB24270 
ORB24280 
ORB24290 
ORB24300 
ORB24310 
ORB24320 
ORB24330 
ORB24340 
ORB24350 
ORB24360 
ORB24370 
ORB24380 
ORB24390 
ORB24400 
ORB24410 
ORB24420 
ORB24430 
ORB24440 
ORB24450 
ORB24460 
ORB24470 
ORB24480 
ORB24490 
ORB24500 
ORB24510 
ORB24520 
ORB24530 
ORB24540 
ORB24550 
ORB24560 


REAL 


INC 


COMMON IWORK( 3800) 


DATA IWDIM/3800/ 

КЕ = 6378105 

SCALE THE EARTH PLOT AND CENTER ON THE ORIGIN 
SCFAC = RE/RSTEP 

SCFAC2 = SCFAC * 2.0 

XPHS = 1525057 5708 US CRAO 

YPHS = 747053205 - СЕ 


ҮРОТЕ = О О ИН ЕЕ г 
ДЕ POLET GI T90) IHREN 
УРОВЕ 2 ТРОШЕ 2290 


ENDIF 


YORIG = YPOLE - 90 


YMAX 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


= YPOLE + 90 
RESET( 3HALL) 

PHYSOR( XPHS, YPHS) 

PROJCT( ' LAMBERT EQ/AREA' ) 
MAPOLE( 0. 0, YPOLE) 

AREA2D (SCFAC2,SCFAC2) 

THKFRM (0. 02) 

GRAF( -90. ,30. ,90. , YORIG,30. , YMAX) 
FRAME 

MAPFIL( 'MAPDTA' ) 
LBLANK('LAND' , IWDIM) 

GRID(1, 1) 

LOLA WATER ТИЕ) 

DASH 

CRID 1.1) 

RESET ('DASH') 

ENDGR(0) 


RETURN 


END 


УУ и У У У је о ел их о ок ede e des 


m 
% 


Р 


SUBROUTINE GRTRK( AINRAY, APRAY , TARAY, TF ,NUM, TIMRAY ) 


DIMENSION AINRAY(500) ,APRAY( 500) , TARAY( 500), 


410 IF (I 
x 
+ 


ще 


ELARAY( 500) , ELORAY( 500) , TLONG( 500) , TLAT( 500) , TIMRAY( 500) 
ВЕ = 6. 3782Е+03 
EROT = 7. 292115856E-05 
STF = (TF) 
I=1 
LOAD ARRAYS WITH LATITUDE AND LONGITUDE 
.LE. NUM) THEN 
- RE*COS(APRAY(I))*COS(TARAY(I))-RE*SIN(APRAY(I))* 
SIN(TARAY(I)) 
Y - RE*COS(AINRAY(I))*SIN(APRAY(I))*COS(TARAY(I)) + 
RE*COS( AINRAY(I))*COS( APRAY( I) )*SIN( TARAY(I)) 
Z 2 RE*SINCAINRAY(I))*SIN(APRAY(I))*COS(TARAY(I)) + 
RE*SIN(AINRAY(CI))*COS(APRAY(I))*SIN(TARAY(I)) 


ВЕ 
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ORB24570 
ORB24580 
ORB24590 
ORB24600 
ORB24610 
ORB24620 
ORB24630 
ORB24640 
ORB24650 
ORB24660 
ORB24670 
ORB24680 
ORB24690 
ORB24700 
ORB24710 
ORB24720 
ORB247 30 
ORB24740 
ORB24750 
ORB24760 
ORB24770 
ORB24780 
ORB24790 
ORB24800 
ORB24810 
ORB24820 
ORB24830 
ORB24840 
ORB24850 
ORB24860 
ORB24870 
ORB24880 
ORB24890 
ORB24900 
ORB24910 
ORB24920 
ORB24930 
ORB24940 
ORB24950 
ORB24960 
ORB24970 
ORB24980 
ORB24990 
ORB25000 
ORB25010 
ORB25020 
ORB25030 
ORB25040 
ORB25050 
ORB25060 
ORB25070 
ORB25080 
ORB25090 
ORB25100 
ORB25110 
ORBZS 20 


430 


CALCULATE LATITUDE 
ELARAY(I) - (ASIN(Z/RE)) * (180/3. 14159) 


TRAP 'X' AND 'Y' FOR ARCTAN IN CALCULATING LONGITUDE 
i COX HE. EOD CANET СЕ ОНО ATHEN 
0 


ELSEIF ((Y .GE. -10). AND. (Y . LE. 0.0)) THEN 


Y = -10. 

ENDIF 

COC БЕТІ МАҚ ТА ТӘН ПӘ) ТИЕ 
Х = 10. 

ЕЗЕН х= СЕ -10). АБ. СХ. БЕ. 0. О) EN 
mecs 1D. 

ENDIF 


CALCULATE LONGITUDE 
ЕОС ИОК А АС АЕ АОЛ TIMRAYCI) 001607 3214159) 


BUDENSIS DES TOS = 06 07 ТОРО) 

TES (EEORAN M L 50) THEN 
ELORAY(I) < ЕТОКАК(Т) + 360 
СОТО 420 

Е ЕЕ 

Iob] 

GOTO 410 

ENDIF 


SET DISSPLA 

БА ШЕК618 

CALL RESET( 3HALL) 

EXER YAXANG СО.) 

CALL PHYSOR(1.0,6.0) 

ШІ МАМЕ 7,1) 

ОДП УКАМЕС 7,1) 

ШЕШІ АҚЕА22(7.5,3.75) 

CALL HEADIN ('GROUND TRACKS‘ ,100,1.5,1) 
GALL, SCMPLX 

CALL MAPGR(-180. ,90. ,180. ,-90. ,30. ,90. ) 
Genk GRID (1,1) 

CALL MAPFIL ('MAPDTA' ) 

m= 


IGNORE Boundary POINTS 
Ше CEGORAY( 1) . LI. LU 
ГЕОС) . СТ. 125). 


(ELARAY(I) .LT. -85) ги 
(ELARAY(I) .GT. 85)) THEN 
ТЕТтАа 1 
GOTO 430 
ENDIF 
ITEMP = 1 


LOAD PlRSt POINT OF NEW PLOT SEGMENT 
ЕЕ ТЕХ 


ORB25130 
ORB25140 
ORE 25150 
ORB25160 
ORB25170 
ORB25180 
ОКЕ ШОО 
ORB25200 
ОКВ 5210 
ORB25220 
ОКВ252 30 
ORB25240 
ORB25250 
ORB25260 
ORB25270 
ORB25280 
ORB25290 
ORB2 35300 
ORB25310 
ORE 5220 
ORB25330 
ОКВ25 340 
BID 
ORB25360 
208255870 
ORB25380 
015256120, 
ORB25400 
ORB25410 
ORB25420 
ORB25430 
ORB25440 
ORB25450 
ORB25460 
ORB25470 
ORB25480 
ORB25490 
ORB25500 
(ИВ ИО 
OREB2Z 5520 
ORB25530 
ORB25540 
ORB25550 
ORB25560 
ORB25570 
ORB25580 
ORB25590 
ORB25600 
ORB25610 
ORB 3670 
ОКВ2 5630 
ОКВ25640 
ORB25650 
ОКВ2 5600 
0КВ2 5670 
ORB25680 


TLONXG(ITEHP) = ELORAI T) 
TLAT(ITEMP) = ELARAY(I) 
ESR ae 
ЕС 1-06. м Е. 
CALL POLY3 
CALL CURVE(TLONG, Tiad , TTEMPS 
ENDIF 
ENDIF 


LOAD SECOND POINT IN LINE SEGMENT 
TFC! «bee NU ne 
Тр el 
TLONGC ITEMP) < ЕТОКАУ(Т) 
TLAT(ITEMP) = ELARAY(T) 


Tm Tp 
ЕНИ О КОН ЕЛ НЕК 
CALL POLY3 


CALL NOCHEK 
CALL CURVE(TLONG,TLAT,ITEMP,1) 
ЕТЕ 
ENDIF 


* LOOP UNTIL SEGMENT REACHES EDGE OR NO MORE POINTS 
440 IF (I .LE. NUM) THEN 


BOTH LAT AND LONG INCREASING 
ТЕГ  CELCKRAMGT ==" 2) SLE Е ТОСКА И ТТ а 


4 (ELARAY(I - 2) .LE. ELARAY(I - 1))) THEN 
IF((ELORAY(1) .LT. -170) .OR. 
+ (ELARAY(I) .LT. -80)) THEN 


CALL POLY3 
CALL NOCHEK 
CALL CURVE(TLONG, TLAT, ITEMP, 1) 
GOTO 430 
ELSE 
TTE - Ро 
TLONGCITEMP) < ЕБОКАУ(Т) 
TEAT(ITEMP JES ELARA T 
ENDIE 


Ы BOTH LAT AND LONG DECREASING 
ELSIF ( (CELORAY(? - 2) GCL ELOR ИНЕСІ 


+ (ELARAY(I - 2) GT. EBARATCI ГЕТЕ 
IFE((ECORAY(I) -GT ПЕ! 
+ (ELARAY( 1) „61, 80)) THEN 
ӨЛЕ аи 


CALL NOCHEK 
CALL СОВУЕСТЕО С, ВАТТ TRR 
GOTO 430 
ELSE 
ТЕМ = ПЕРА 
TLONG(ITEMP) = ЕГОКАУ(1) 
LAT(ITEMP) = ЕГАКАХ (1) 
Е 2ТЕ 


LAT INCREASING, LONG. DECREASING 


ORB25690 
ORB25700 
ORB25710 
ORB25720 
ORB25730 
ORB25740 
ORB25750 
ORB25760 
ORB25770 
ORB25780 
ORB25790 
ORB25800 
ORB25810 
ORB25820 
ORB25830 
ORB25840 
ORB25850 
ORB25860 
ORB25870 
ORB25880 
ORB25890 
ORB25900 
ORB25910 
ORB25920 
ORB25930 
ORB25940 
ORB25950 
ORB25960 
ORB25970 
ORB25980 
ORB25990 
ORB26000 
ORB26010 
ORB26020 
ORB26030 
ORB26040 
ORB26050 
ORB26060 
ORB26070 
ORB26080 
ORB26090 
ORB26100 
ORB26110 
ORB26120 
ORB26130 
ORB26140 
ORB26150 
ORB26160 
ORB26170 
ORB26180 
ORB26190 
ORB26200 
ORB26210 
ORB26220 
0882620 
ORB26240 


elo e o mao n c m'a 


ВЕБЕ ВОВА (ГГ жазу СТ ЕВ му) 


TAN. 


(ELARAY(I - 2) .LE. ELARAY(I - 1))) THEN 
IF((ELORAYCI) .GT. 170) .OR. 
(ELARAYCI) .LT. -80)) TEEN 
CALL POLY3 
CALL NOCHEK 
CALL CURVE(TLONG ,TLAT, ITEMP, 1) 
GOTO 430 
ELSE 
[TEMP = ITEMP + 1 
TLONG(ITEMP) = ELORAY(I) 
TLAT( ITEMP) = ELARAY(I) 
ENDIF 
LAT. DECREASING, LONG. INCREASING 
ELSEIF((ELORAY(I - 2) .LE. ELORAY(I - 1)) . AND. 
(ELARAY(I - 2) .GT. ELARAY(I - 1))) THEN 
IF((ELORAY(I) .LT. -170) .OR. 
(ELARAY(I) .GT. 80)) THEN 
CALL POLY3 
CALL NOCHEK 
CALL CURVE(TLONG,TLAT, ITEMP, 1) 
GOTO 430 
ELSE 
ITEMP = ITEMP + 1 
TLONG(ITEMP) = ELORAY(I) 
TLAT(ITEMP) = ЕГАВАУ( Г) 
ENDIF 
ENDIF 
IF( I .EQ. NUM) THEN 
CALL POLY3 
CALL NOCHEK 
CALL CURVE(TLONG,TLAT, ITEMP, 1) 
ENDIF 
mo 
GOTO 440 
ENDIF 
CALL POLY3 
CALL NOCHEK 
CALL CURVE(TLONG,TLAT, ITEMP, 1) 
CALL COMPLX 
CALL HEIGHT(.2) 
CALL THKFRM (0. 03) 
CALL FRAME 
CALL RESET('COMPLX') 
CALL RESET( 'HEIGHT') 
CALL ENDGR (0) 
RETURN 
END 


рт реа 


. . з . 
l'on ann n on e n un on a n. 


ORB26250 
ORB26260 
ORB26270 
ORB26280 
ORB26290 
ORB26300 
ORB26310 
ORB26320 
(Бов о О 
ORB26340 
ORB26350 
ORB26360 
ORB26370 
ORB26380 
ORB26390 
ORB26400 
ORB26410 
ORB26420 
ORB26430 
ORB26440 
ORB26450 
ORB26460 
ORB26470 
ORB26480 
ORB26490 
ORB26500 
ORB26510 
ORB26520 
ORB26530 
ORB26540 
ORB26550 
ORB26560 
ORB26570 
ORB26580 
ORB26590 
ORB26600 
ORB26610 
ORB26620 
ORB26630 
ORB26640 
ORB26650 
ORB26660 
ORB26670 
ORB26680 
ORB26690 
ORB26700 
ORB26710 
ORB26720 
ORB26730 
ORB26740 
ORB26750 
ORB26760 
ORB26770 
ORB26780 
ORB26790 
ORB26800 


SUBROUTINE DATA( I.A, E,TFEA,TFSU, TF0, TFDRA PER PI TOTT LAR O ORB26810 


+ TDMM,TDMA, TDLAN,TDH, TDAP,MM,MA , LAN,H, AP, V, R) ORB26820 
THIS SUBROUTINE Displays THE ORBITAL DATA FOR BOTH THE PERIFOCAL ORB26830 
AND THE GROUND TRACK PLOTS. ORB26840 
REFER TO DISSPLA USER'S MANUAL FOR EXPLANATION OF DISSPLA ORB26850 
SUBROUTINES ORB26860 

ORB26870 
REAL I,MM,MA, LAN ORB26880 

ORB26890 
MU = 3. 986012E+05 ORB26900 

ORB26910 
CALCULATE THE AVERAGE FORCES FROM THE TOTAL MAGNITUDE OF ORB26920 
FORCE CHANGES ORB26930 
AVGFE = TFEA/50. 0 ORB26940 
AVGFS = TFSU / 50.0 ORB26950 
AVGFM = TFMO / 50.0 ORB26960 
AVGFD = TFDRA / 50.0 ORB26970 

ORB26980 
CALCULATE ORBITAL ELEMENTS IN Usable UNITS ORB26990 
PERH = PER/3600 ORB27000 

ORB27010 
DI = I * (180.0/PI) ORB27020 
DLAN = LAN * (180.0/PI) ORB27030 
DAP = AP * (180.0/PI) ORB27040 

ORB27050 
CALCULATE Average CHANGE IN ELEMENTS FOR ONE PERIOD ORB27060 
АВЕО 50S ОЕВ27070 
AVGDA = TDA / 50.0 ORB27080 
AVGDE = TDES SNO ORB27090 
AVGDMM = TDMM / 50.0 ORB27100 
AVGDMA = TDMA / 50.0 ORB27110 
AVGLAN = TDLAN / 50.0 ORB27120 
AVGDH = ТОН / 50.0 08827130 
AVGDAP = TDAP / 50.0 ORB27140 

ORB27150 
CALCULATE RADIUS. SANDE BLOG hina 5 ORB27160 
ЕМЕ 2 ((V**2)/2) - (MU/R) ORB27170 
ВР = А*(1 - E) 08827180 
RA = A*(1 + E) ORB27190 
VP - SQRT(2*(ENR + (MU/RP))) ORB27200 
VA = SQRT(2*(ENR * (MU/RA))) ORB27210 

ORB27220 

ORB27230 
SET DISSPLA ORB27240 
CALL RESET( 3HALL) ORB27250 
CALL CUPL ORB27260 
CALL PHYSOR(COSONONMGD ORB27270 
CALL AREA2D(9 SA0) ORB27280 

ORB27290 
PRINT DATA ORB27300 
CALL MESSAG('I s 5" ,100,0. 25,3. 67) ORB27310 
CALL REALNO(DI 3, ABUT ТАРЫ ORB27320 
CALL MESSAG(' DEG.$',100, ABUT , ABUT') ORB27330 
CALL MESSAG( А- $1100, ABUT PADCI) ORB27340 
CALDMREALNOCA, 1. АВИТ а ORB27350 
CALL MESSAG(’ KMS$',100, ABUT’ , АВИТ") ORB27360 


“4 


CALL 
CALL 
CALL 
CALL 
CALL 


CALL 


CALL 


MESSAG(' E = $S T100, ABUT ' ABUT' ) 
REALNO(E,3, ABUT’ , 'ABUT’) 

MESSAG(' PER = $',100,'ABUT' ,'ABUT') 
REALNO( PERH,2,'ABUT' ,'ABUT' ) 

MESSAG(' HOURSS' ,100, 'ABUT' ,'ABUT') 


MESSAG(.' AVERAGE RATE OF CHANGE OF ELEMENTS PER SECOND $', 
100,1.0,3.0) 


MESSAG('DI/DT = $',100,0. 25,2. 67) 
REALNO(AVGDI, -2, ' ABUT' , ' ABUT!) 

МЕЗО DA/DT = $100 "ABUT , ABUT ) 
REALNO( AVGDA,-2,'ABUT','ABUT') 

MESSAGO  —DE/DT s $ OO ABUT | АВОТ ) 
REALNO( AVGDE ,-2,' ABUT’ , 'ABUT' ) 


MESSAG( 'DMM/DT = $' ,100,0. 25,2. 33) 
REALNO( AVGDMM,-2, ABUT','ABUT') 

MESSAG(' DMA/DT 7 $',100,'ABUT', ABUT ) 
REALNO(AVGDMA, -2, ABUT' , ABUT' ) 

MESSAG(' DLAN/DT = $',100,'ABUT' ,' ABUT’ ) 
REALNO(AVGLAN,-2, ABUT' , ABUT') 


МЕББАСГ ТЕН ТЕ Е- 52100 0225 200) 
REALNO( AVGDH,-2,'ABUT' ,'ABUT') 

MESSAG(' DAP/DT = $',100,'ABUT' ,'ABUT') 
REALNO( AVGDMA , -2, ' ABUT' , ' ABUT' ) 


MESSAG('AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KM/S**2) 


2:0 100. 1.0, 1. 67) 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


CALL 
CALL 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


MESSAGO EARTH-- $ 100,0. 10,1. 33) 
REALNO(AVGFE,-1, ABUT' , АВОТ") 
MESSAG(' MOON 2» $',100,'ABUT' , ' ABUT' ) 
REALNO(AVGFM,-1,'ABUT', ABUT') 
MESS46, 505 = 5.100. ABUT, ABUT’) 
REALNO(AVGFS,-1,'ABUT' , ' ABUT' ) 
MESSACOMDRAG S S 100. ABUT. - ABUT ) 
REALNO( AVGFD,-1,'ABUT' , ABUT') 


MESSAG( ' PERIGEES' ,100,2. 75,1. 0) 
MESSAG( ' Apogee$',100,' ABUT' , ABUT ) 


MESSAG('RADIUS (KM)$',100,0.25,0.67) 
MESSAG('RP -2$',100,2. 75,0. 67) 
REALNO(RP,1,'ABUT','ABUT') 

MESSAG(' © ОТОО, ABUT: ABUT ) 
MESSAG(' RA =$',100,'ABUT’, ABUT ) 
REALNO(RA,1,'ABUT' ,' ABUT' ) 


HESSAC VELOCITY (KM/SECG)S , 1000. 25,0. 33) 
MESSAG('VP 2$',100,2. 75,0. 33) 
REALNO(VP,2, ' ABUT' ,' ABUT' ) 

MESSAG(' СОО ДВЕ. АВО) 
MESSAG(' VA -$',100,'ABUT' ,' ABUT ) 

ВВ И.о. АВ, ABUT ) 


~ 
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ШЕВ 370 
ORB27380 
ORB27390 
ORB27400 
ORB27410 
ORB27420 
ORB27430 
ORB27440 
ORB27450 
ORB27460 
ORB27470 
ORB27480 
ORB27490 
ORB27500 
ORB27510 
ORB27520 
ORB27530 
ORB27540 
ШЕВО 50 
ORB27560 
ОКВ2 7570 
ОКВ2 7580 
ORB27590 
ORB27600 
ORB27610 
ORB27620 
ORB27630 
ORB27640 
ORB27650 
ORB27660 
ORB27670 
ORB27680 
ORB27690 
ORB27700 
ORB27710 
ORB27720 
ORB27730 
ORB27740 
ORB27750 
ORB27760 
ORB27770 
ORB27780 
ORB27790 
ORB27800 
ORB278 10 
ORB27820 
ORB27830 
ORB27840 
ORB27850 
ORB27860 
ORB27870 
ORB27880 
ORB27890 
ORB27900 
ORB27910 
ÜRD27920 


CALL RESET('COMPLX') 
CALL ENDGR(0) 

RETURN 

END 


ОКВ? 
ORB27 
ORB27 
ORB27 
ORB27 
ORB27 


APPENDIX B. COORDINATE SYSTEMS 


А. ‘IJK’: GEOCENTRIC - EQUATORIAL 


a 


vernal equinox 
direction 





Figure 3. Geocentric-equatorial coordinate system 


The geocentric-equatorial system as seen in Figure 3 has its origin at the earth’s 
center. The fundamental plane is in the equator and the positive X-axis points in the 
vernal equinox direction. The Z-axis points in the direction of the north pole. This 
svstem 1s not fixed to the earth and turning with it; rather, the geocentric-equatorial 
frame 1s nonrotating with respect to the stars (except for precession of the equinoxes) 
and the earth turns relative to it. Unit vectors. /. J, and K shown in Figure 3. lie along 
m Y. and Z respectively. [Ref. I: p.55] 


~, 
~, 


В. ‘’PQW’: PERIFOCAL 





Figure 4. X Perifocal coordinate system 


The perifocal coordinate system has its fundamental plane in the plane of the satel- 
lite’s orbit as seen in Figure 4. The coordinate axis are named, A.. Y, and Zee 
X, axis points toward the perigee; the Y, axis 15 rotated 90 degrees in the direction of 
orbital motion and lies in the orbital plane; the Z, axis along h completes the right- 
handed perifocal system. Unit vectors in the direction of X,, Y,, and Z, are called P, Q 


апа respectively Ке Ев | 
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С. ‘RSW’: ORBITAL 





Figure 5. Orbital coordinate system 


(Figure 9.4-1, Ref. 1) 

The orbital coordinate system has its principle axis, R (unit vector r), along the in- 
stantaneous radius vector, r as seen in Figure 5. The axis S is rotated 90 degrees from 
К in the direction of increasing true anomaly. The third axis. W, is perpendicular to 
both R and S. Note that this coordinate system is simplv rotated v, from the PQW 
perifocal svstem. [Ref. 1: p.39$] 


D. COORDINATE TRANSFORMATIONS 

The coordinate transformations, for the previous coordinate svstems. use angular 
rotations about the axis to evaluate the transformation matrix. The matrix elements r, 
are calculated, then applied to the old vector to get the vector in the new coordinate 
system. The following orbital elements are used: 


Q = longitude of ascending node 


со = argument of perigee 
1 = inclination 

ш = argument of latitude 
Vo = true anomaly 


The coordinate transformations follow [Ref. 1: p.74-83] 


РОМ СОМИ 
ra= cos coso =~ sin S sin OKOS 
а= = СОЗ ОИ о = е O Ges cos 


оо ао о 


г = п О сох E cos SIEG OS 
r= — $in Q sin w + cos Q cos œ cos i 
y4 — — cos Q sini 

[seen ca ST 

ка = COS Q SIMI 


ju esr 
ЕД” 
Пп OPE 
К 


--- 


ва? + 0 = T 


ато Ооа 
Е ди + zu + ғы 
O = ns + rad + rk 
= = E 2 + К 
Sy WK tomo 


y "cos GOS 1 = SNU COS 
го = ЗО со ио SIDES COM. 


ЖЕСІ СПП 
Fa = — COS Q SIN ty - Sin Q cos u, cos i 
= Г О) ып ну COS © COS mace 


СТЕ! 
а 


pc cO n 

"у = COS! 

R= А ES 
2. ПИ 
H E + А Lm 


КИ то LJ Ko (inverse of =3) 
í = mo -- S - ET ІШ 

J = rR S cS 

и 


~ 


oe 


ао оа T + 3 
г 
meee OW to RSW 
Ку = COS Vo 
=з, 


г; = 0.0 

И 

собе 

Eu 

ка + 0.0 

3; = 0.0 

ка = 1.0 

К =>"; y — О - пз 
У = ИР в. + ий” 
а РО + 


6. RSW o PQW (inverse of 25) 
= ғ ns == EN 
tnt о: 
г. + до i nV 
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АРРЕХФТА С. ОКВТТАГ ЕГЕ УК 5 


The user is assumed to be studving orbital mechanics and should understand the 


orbital elements and how to calculate them. A brief description of the elements and the 


equations used to calculate the elements follow. For a detailed explanation of the ele- 


ments and the equations to calculate them refer to Chapters 1 and 2 of reference 1. 


Figure 6 on page 83 shows the orbital elements in. the Geocentric-Equatorial and 


perifocal coordinate svstem. 


i 


„ә 


ӨТ 


Angular Momentum (В): 
The specific angular momentum is a constant of the motion of the satellite. 


Че те а ле 
роже = АГ НАУ + А, К 


p 


= к= 


› 1, 
/ ғ 


hj — ry; — riy 


Е 
& 3 / t 


- 


h= dp rhs hy 


Nodes secon iim: 
The node vector is a vector poinung along the line of nodes in the direction 
of the ascending node. 


— amb 


п = Кх = = + hj 
n=, h +h 


Semi-latus rectum (p): 
The semi-latus rectum is a geometric constant of the conic section. 


3 
А“ 
í 


fct um 


Ессепішете е) 
The eccentricity is a constant defining the shape of the conic orbit. 
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Figure 6. Orbital elements 


The semi-major axis is a constant defining the size of the orbit. 


i-e 


а - р 


6. Inclination (1): 


The inclination is the angle between the ‘K’ unit vector in the ‘IJK’ system and 
the angular momentum Vector, ‘h’. 
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~ 
. 


19. 


ide 


— ым 


я —:, ГА) А —] hy, 
jo ML = о о 

{ fl 
Longitude of ascending node (22): 

The longitude of the ascending node 1s the angle in the fundamental plane , 
between the ‘]’ unit vector and the point where the satellite crosses throughgmubs 
fundamental plane in a northerly direction (ascending node) measured counter- 
clockwise when viewed from the north side of the fundamental plane. 


( --1 п; 
MOS. 
Argument of perigee (@): 

The argument of perigee is the angle in the plane of the satellite’s orbit, be- 
tween the ascending node and the perigee point, measured in the direction of the 
satellite's moton. 

па о) 


oO = cos ( = = 


ne не 


True anomaly at epoch (v,): 

The true anomaly at epoch is the angle in the plane of the satellite's orbit. be- 
tween perigee and the position of the satellite at a particular time. 4. called the 
“epoch’. 


ae саћа 
^ . 


=l; €! 
СО т | 
Argument of latitude (u,): 
The argument of latitude is the angle in the plane of the orbit. between the 
ascending node and the radius vector to the satellite at time 1. 


о —_—а 


oly НР 
He COS) Ga) 
Longitude of perigee (I1): 

The longitude of perigee is the angle from ‘Il’ to perigee measured eastward to 
the ascending node and then in the orbital plane to perigee. 


П-0-- о 


12. True longitude at epoch (4): 


The true longitude at epoch is the angle between ‘1’ апа ғ, (the radius vector 
to the satellite at 4, measured eastward to the ascending node and then in the orbital 
plane to 7. 


152 РеПОШ ЕТ) 


14. 


The period is the time the for the satellite to complete one orbit. 


ха 


Бер = В г! 


а EAN: 
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Шише сесешиле anomal is the вое Cetween ithe percee and a position on ani 
auxiliary circle circumscribed about the ellipse where a perpendicular line to the 
major axis has been extended from the epoch location of the satellite to the auxul- 
E circle. 

е + cos(v) 


К = cos | RE те 
| + e costr) 


Ealan motion (n): 
The mean motion is defined below: 





М’ = [= 
“а? 


16. Меап апотају (МА): 
The mean anomialv 1s defined below: 


МА =п( = 1) = ЕЛ — е зи(ЕА) 


Eime of flight (TF): 
The time of flight is the elapsed time from when the satellite was at perigee to 
Mae current epoch. 
г 


га 


Nl 1) = \/ “ДЕ (EA RC sin( EA )) 
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APPENDIX D. SAMPLE ORBITS 


To demonstrate the capabilities of the program, a variety of orbital plots will follow: 
l. Low earth orbit (LEO) 


PERIFOCAL COORDINATE SYSTEM 


1 = 45.000 ОРС. А = 7222.2 КМ Е = 0.100 РЕН = 1.70 HOURS 


AVERAGE RATE OF CHANGE OF ELEMENTS PER SECOND 
DIYDT = 0.00 DA/OT = 0.00 DE/DT = 0.00 
DMM/DT = 0.00 DMA/DT = 0.00 DLAN/DT = 0.00 
DH/DT = 0.00 DAP/DT = 0.00 

AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KM/S°°2) 

EARTH = 0.0 MOON = 0.0 SUN= 0.0 DRAG = 0.0 
PERIGEE Apogee 

RADIUS (KM) RP = 6500.0 ВА = 7944.4 
VELOCITY (KM/SEC) УР = 821 VA = 6.72 





Figure 7, Unperturbed Low Earth Orbit (LEO) 


Figure 7 shows the perifocal plot of a satellite in an unperturbed low earth 
orbit (LEO). The initial parameters of the orbit were entered as follows: 
radius of perigee (RP) = 6500 km 
eccentricity (e) = 0.1 
inclination (1) = 45 degrees. 
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PERIPOCAL COORDINATE SYSTEM 


1 = 44.998 DEG. A= 7203.3 КМ Е = 0.088 PER= 1.69 HOURS 


AVERAGE RATE OF CHANGE OP ELEMENTS PER SECOND 
DI/DT = 4.20°10°° DA/DT = 8.36°10°" DE/DT = 1.007107 
DMM/DT = 1.80°10°° DMA/DT = 9.21°10`° DLAN/DT = 9.48'10"" 
DH/DT = 5.83°10°° DAP/DT = 9.21°10°° 
AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KM/S*72) 
EARTH = 8.6*10'* MOON = 9410“ SUN » 4.3"10'" DRAG = 1.4°10`* 
PERIGEE Apogee 
RADIUS (KM) RP = 6499.6 RA = 7907.0 
VELOCITY (KM/SEC) УР = 8.20 VA = 6.74 





Figure 8. Perturbed Low Earth Orbit (LEO) 


With perturbing forces applied to the previous LEO, the drag force will be the 
dominate perturbing force. The drag will act as a negative velocitv change applied 
in the area of perigee, with the result of decreasing the semi-major axis length, this 
in effect will decrease the eccentricity of the orbit, as can be seen bv comparing the 
orbital data of the unperturbed LEO in Figure 7 on page 86 with the orbital data 
phe perturbed LEO in Figure 8. 
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2. Circuli Orbit. 


GROUND TRACK 


| = 59.967 DEG. A= 6986.1 KM E= 0.000 PER= 1.61 HOURS 


AVERAGE RATE OF CHANGE OF ELEMENTS PER SECOND 
DI/DT = 4.02'10' DA/DT - 9.2410 РЕ/ОТ = 5.93"10"" 
рмм/рт = 2.26710 * DMA/DT = 1.52*10'' DLAN/DT + 7.32"10” 
DH/DT = 5.86'10' DAP/DT = 1.52710" 
AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KM/S*°2) 
EARTH = 1.2°10°° MOON = 8.8°10°'° SUN = 4.3°10°'° DRAG = 3.0°10 
PERIGEE Apogee 
RADIUS (KM) RP = 6984.4 RA = 6987.8 
VELOCITY (KM/SEC) VP = 7.56 VA = 7.55 





Figure 9. Circular Orbit 


An example of the plot of the ground track of a sequence of three 60 degree 
inclined perturbed circular orbits with a radius of 7000 km is shown in Figure 9. 
The sequence of orbits displays the precession of the orbit around the earth. 
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Beri tcansier orbit. 


PERIPOCAL COORDINATE SYSTEM 


220216 


220215 
= 0.000 DEG. A= 14510.8 KM E= 0.518 PER= 4.83 HOURS 


AV ERAGE RATE OF CHANGE OF ELEMENTS PER SECOND 
DI/DT = 0.00 DA/DT = 0.00 DE/DT = 0.00 
DMM/DT = 0.00 DMA/DT = 0.00 DLAN/DT = 0.00 
DH/DT = 0.00 DAP/DT = 0.00 
AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KM/S**2) 
EARTH = 0.0 MOON = 0.0 SUN = 0.0 DRAG = 0.0 


PERIGEE Apogee 
RADIUS (KM) RP = 7000.0 RA = 22021.5 
VELOCITY (KM/SEC) УР = 9.30 УА = 2.95 





Figure 10. Transfer Orbit 


The transfer orbit between a circular, equatorial LEO and a molniva orbit 
(high eccentric orbit) is shown in Figure 10. A velocity increase of 1.75 km's was 
applied at the perigee to simulate a perigee kick to boost the satellite into the 
molniva orbit. A similar velocity change could then be applied at apogee to create 
a high altitude circular orbit. or a negative velocity change applied at perigee could 
be used to bring the satellite back to a LEO. 
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4. Geosvnchronous orbit 


GROUND TRRCK 


60.000 DEG. A 7 42234.5 KM E- 0.000 PER - 23.99 HOURS 


AVERAGE RATE OF CHANGE OF ELEMENTS PER SECOND 
DI/DT - 1.0410 * ОА/ОТ = 1.26“10“ DE/DT = 1.08°10°° 
DMM/DT - 326"10'' DMA/DT - 3.02*10* DLAN/DT = 1.89710" 
ОН/ОТ = 3.55*10“ ОАР/ОТ = 3.02“10“ 
AVERAGE MAGNITUDE OF FORCES PER UNIT MASS (KN/S°°2) 
EARTH = 9.3°10°° MOON = 4.9°10°° SUN = 2.6°10°° DRAG = 3.110“ 
PERIGEE Apogee 
RADIUS (KM) RP = 42233.9 RA = 42235.3 
VELOCITY (KM/SEC) VP = 3.07 VA = 3.07 





Figure 11. Geosynchronous Orbit 


The ground track of a perturbed geosynchronous orbit inclined 60 degrees is 
shown in Figure 1l. The orbit displays the figure eight typical with inclined 
geosvnchronous orbits. 
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